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ABSTRACT 


This  thesis  presents  a  general  development  of  continuous  piece- 
wise  interpolation  formulae  extended  in  two  directions;  first  the 
extension  requiring  continuity  of  higher  order  derivatives;  and 
secondly,  the  extension  in  the  multi  step  direction.  An  algorithm 
and  a  general  program  for  the  computation  of  the  coefficient  matrices 
of  the  continuous  piecewise  interpolation  formulae  in  floating-point 
form  are  also  given.  These  formulae  are  tested  and  compared  on 


known  functions. 
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CHAPTER  I 


INTRODUCTION 

Interpolation  has  been  described  as  the  art  of  reading  between 
the  lines  of  a  table.  In  elementary  mathematics  the  term  usually 
refers  to  the  process  of  computing  intermediate  values  of  a  function, 
such  as  trigonometric  and  logarithmic,  from  a  set  of  tabular  values 
of  that  function.  In  engineering  the  function  is  usually  unknown  but 
is  given  in  terms  of  a  table  of  values  obtained,  for  example,  from 
measurements  of  physical  quantities  such  as  stream  flow,  precipita¬ 
tion,  etc..  For  convenience  the  measurements  may  be  made  at  equal 
intervals  of  an  independent  variable.  According  to  Weierstrass 's 
theorem  the  unknown  function  can  be  approximated  by  a  polynomial  to 
any  desired  degree  of  accuracy  provided  that  it  is  continuous  over 
a  finite  interval  of  interest  and  sufficient  data  is  available.  How¬ 
ever,  in  some  engineering  problems  the  use  of  a  single  polynomial  for 
the  whole  set  of  data  points  may  be  impractical.  One  reason  is  that 
when  the  data  set  is  very  large  the  interpolation  process  may  be  too 
time-consuming  if  not  also  inaccurate  due  to  round-off  error.  Another 
reason  is  that  the  data  set  may  be  open-ended,  e.g.,  data  is  collected 
continually  every  day.  Therefore  in  some  engineering  problems  it  is 
desirable  or  even  necessary  that  interpolation  be  applied  to  a  subset 
of  data  points  at  a  time.  The  process  may  thus  be  called  "piecewise 
interpolation". 


. 
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The  derivative  of  the  function  in  engineering  problems  usually 
represents  a  physical  quantity.  Therefore,  it  is  undesirable  to  have 
two  values  of  the  derivative  at  the  joining  point  between  any  two 
interpolating  polynomials.  The  process  of  piecewise  interpolation 
with  continuous  derivatives  may  thus  be  called  "continuous  piecewise 
interpolation". 

1 . 1  History 

The  theorem  of  Weierstrass  was  first  proven  in  1885  and  since 
then  a  vast  amount  of  literature  on  classical  interpolation  has  been 
written.  Excellent  discussions  of  the  basic  concepts  of  interpola¬ 
tion  may  be  found  in  Hilderbrand  [17],  Kopal  [21]  and  Scarborough  [25]. 
The  discussions  on  classical  interpolation  as  found  in  these  and  other 
standard  texts  on  numerical  analysis  deal  with  concepts  of  "piecewise 
interpolation"  without  any  reference  to  the  process  of  "continuous 
piecewise  interpolation".  An  interpolation  process  closely  related 
to  continuous  piecewise  interpolation  is  the  well  known  spline 
interpolation. 

The  theory  of  spline  approximation  first  appeared  in  a  paper  by 
Schoenberg  [26].  A  spline  function  of  degree  (m  -  1),  having  the  knots 
at  x  =  1,2,...,  (n  -  1),  is  a  function  of  the  class  °°) 

which  reduces  to  a  polynomial  of  degree  at  most  (m  -  1)  in  each  of  the 
intervals  (-°°,  1 ) ,  (1 ,  2) ,  . . . ,  (n  -  1 ,  «)•  A  monospline  K(x),  of 
degree  m,  may  be  obtained  if  we  add  the  monomial  xm/m!  to  a  spline 


. 


3 


function.  Much  of  the  present-day  theory  of  splines  as  based  on 
Holloday's  Theorem  and  its  proof  [18],  requires  that  the  second  deriva¬ 
tive  at  the  junction  point  between  the  two  interpolation  polynomials 
be  equal  to  zero.  This  will  not  guarantee  that  the  function  will  be 
continuous  at  the  junction  point. 

After  1946,  Schoenberg,  together  with  some  of  his  students,  con¬ 
tinued  investigations  of  splines  and  monosplines.  In  particular, 
Schoenberg  and  Whitney?[27]  and  [28],first  obtained  criteria  for  the 
existence  of  certain  splines  of  interpolation.  For  the  case  of  splines 
of  even  order  with  interpolation  at  the  junction  points  a  simpler 
approach  to  the  question  of  existence  was  given  by  Ahlberg,  Nilson, 
and  Walsh_,[4]  and  [6].  The  convergence  of  the  spline  approximations 
has  also  come  under  close  scrutiny.  A  detailed  analysis  was  presented 
by  Ahlberg  and  Nilson  in  [1],  [2]  and  [3],  with  a  modification  of 
the  restriction  by  Sharma  and  Meir  in  [29]  and  [30].  The  theory  of 
splines  has  also  been  extended  in  a  number  of  directions,  in  particular 
the  extension  to  several  dimensions  by  de  Boor3[12]  and  [13],  and  later 
by  Ahlberg,  Nilson,  and  Walsh  in  [5]  and  [7], 

Even  though  many  papers  have  been  written,  since  the  introduction 
of  spline  approximation  theory  in  1946,  no  reference  is  found  to  the 
process  of  "continuous  piecewise  interpolation"  as  introduced  in  1961 
by  Snyder  [32]  for  a  polynomial  of  degree  up  to  two  with  continuous 
first  order  derivatives;  and  extended  by  the  same  author  in  1967  [34] 
to  polynomials  of  degree  up  to  four  with  continuous  first  and  second 
order  derivatives. 
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1 .2  Purpose  of  the  Study 

Both  of  Snyder's  formulae  may  be  classified  as  "single-step" 
since  the  interpolating  polynomial  is  moved  to  the  right  one  data  point 
at  a  time.  In  this  paper  we  present  a  general  development  for  formulae 
extended  in  two  directions:  firstly,  an  extension  to  higher  order  de¬ 
rivatives;  and  secondly,  an  extension  in  the  multistep  direction.  Also 
presented  is  a  general  computer  program  for  the  computation  of  the  coef¬ 
ficients  of  the  matrix  equation  defining  the  various  interpolation 
formulae. 


. n  •  ...  ro i  ]  upa  i  9  ■}  :•»  a:hu-i.i 


CHAPTER  II 


BASIC  THEORY 


2.1  Polynomial  Interpolation 

Let  f (x )  be  a  function  which  may  or  may  not  be  known  in  closed  form 
but  the  values  of  f^(x.j  )^  where  the  superscript  j  indicates  the  j*'*1 
derivative,  are  given  at  the  (n  +  1)  points,  i.e. 

(2.1)  f^(Xj)  =  yi j  i  =  0,1,... ,  n 

j  =  0,1,...,  (ai  -  1) 

for  some  integers  a.  and  reals  y...  The  problem  of  interpolation  is  to 

I  I  J 

approximate  f(x)  by  another  function  <J>(x)  constructed  such  that 

(2.2)  <f>^(x.)  =  y..  i  =  0,1,...,  n 

j  ~  0,1,..,,  (a..  -  1 ) 

The  interpolating  function  <f>(x)  may  take  a  variety  of  forms  such  as  a 
polynomial,  an  exponential  function,  a  trigonometric  function,  a  con¬ 
tinued-fraction  expansion,  etc..  Since  a  polynomial  can  be  easily 
evaluated  on  a  digital  computer  it  will  be  the  only  form  used  in  this 
thesis.  A  polynomial  of  degree  n  is  defined  as 

Pn<X)  =  X  9i  ^ 


(2.3) 


\  ■  •  ‘  .  • 
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where  n  is  a  positive  integer  and  a.  for  i  =  0,1,,,.,  n  are  constants 

such  that  a„  is  non-zero, 
n 

2.1.1  Notation  and  Basic  Definition  Upper-case  Latin  letters 
designate  matrices  while  all  lower-case  letters  specify  scalars.  A 
rectangular  matrix  A  =  (a..)  with  m  rows  and  n  columns  is  said  to  be 

'  J 

of  dimension  m-by-n.  If  m  =  n,  the  matrix  A  is  square  and  of  order  n. 
A"^  designates  the  inverse  of  a  nonsingular  matrix  A  and  |A|  the  det¬ 
erminant  of  A.  Derivatives  are  specified  by  a  superscript  in  paren¬ 
thesis,  i.e.,  f^(x)  is  the  derivative  of  f(x).  The  largest 
integer  <  £  will  be  denoted  by  [£]  . 

2.1.2  The  Weiers trass  Approximation  Theorem  The  justification 
for  replacing  a  given  function  by  a  polynomial  representation  rests 
on  the  theorem  stated  by  Weiers trass  in  1885. 

Theorem  2.1  lb  ^(x)  Is  continuous  on  a  {^Invtc  interval  [a,b], 
then,  given  any  e  >  0,  thene  extsts  an  n  =  n(e)  and  a  polynomial 
P  ( x )  ob  degnee  n  suck  that  | ^ (x)  -  P  (x)  |  <  e  bon  all  x  In  [a, 6]. 

The  proof  by  Bernstein  of  this  theorem  may  be  found  in  Ralston  [24]. 

This  theorem  justifies  the  use  of  polynomial  approximations  since  it 
guarantees  that  a  polynomial  can  be  found  with  an  arbitrary  small  maxi¬ 
mum  deviation  from  f(x)  on  [a,b]. 


t: 

■ 

: 
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2.2  Hermite  Interpolation  Polynomial 

Suppose  there  is  a  basic  system  of  interpolation  functions 
4>0(x),  4>-j  (x) ,  ...  ,  *  (x),  ...  on  [a,  b].  It  is  required  to  find 
a  linear  combination 


m 

(2.4)  <|>(x)  =  l  c.  <j>. (x)  , 

i=0  1  1 

n 

such  that  (2.2)  is  satisfied.  Since  Y  a.  conditions  have  been  im- 

i=0  1 

posed  on  4>(x )  the  problem  will  always  have  a  unique  solution  if 


(2.5) 


and 


(2.6) 


m  =  (  l  a.)  -  1 
i=0  1 

♦l<xo) .  *m(xo) 

*;<x0>  *i(xo) . ^(xo) 

(a  - 1)  (a  -1)  (an“ 1 ) 

*o  (xo}  *1  (xo}  .  ^m  (xo> 

^0(xi  ^  W . w 

(a-1)  (a-1)  (a  -1 ) 

^o  ^xn^  *1  ^xn^  .  ^m 


t  0 


.  .  '  ; 

•  .« ( * )  f  4 

■ 


. 
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Thus  it  is  possible  to  construct  an  algebraic  polynomial  Hm(x)  of  de¬ 
gree  m  or  less  which  satisfies  the  specified  conditions.  This  alge¬ 
braic  polynomial  may  be  referred  to  as  the  Hermite  interpolation 
polynomial  [9] . 

2.2.1  Uniqueness  of  the  Hermitian  System  We  first  introduce  the 
Alternative  Theorem  [11],  the  application  of  which  is  sufficient  to  show 
that  an  interpolation  problem  has  a  solution,  then  we  apply  this  theorem 
to  the  Hermitian  system  to  show  uniqueness. 

n 

Theorem  2.2  The  homogeneous  system  l  a..x.  -  0 

5=0  ^  i 

l  =  0,1,...,  n  possesses  a  non-txlvlal  solution  U.e.  a  solution  othex 

than  x .  -  0  £ on.  all  l  =  0,1,...,  n)  l{  and  only  ifi  |A|  =  0.  \ ^ox  a 

given  A  =  [a.  .)  thexe  axe  solutions  to  the  non- homogeneous  system,  l.e., 
n  ^ 

l  a-  ■  x  -  =  b .  I  =  0, 1 ,...,  n,  fiox  evexy  selection  ofi  the  quantities 
j-Q  ^5  5  A- 

b.,  then  |A|  i  0  and  the  homogeneous  system  has  only  the  txlvlal  solution. 

We  now  show  that  if  P(x)  is  a  polynomial  of  degree  N  and  satisfies 

(2.7)  P(J)(xi)  =0  for  1  =  0,1,...,  n 

and  j  =  0,1 . (m.-l) 

n 

where  N  =  (  7  m.)  -  1,  then  P (x )  must  vanish  identically.  By  the 
i=0  1 

Factorization  Theorem  [11], 


1  f  ? 

" 

> 


•  . 

. 


, 
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Theorem  2.3  F  (x)  u>  a  polynomial.  o{,  degree  n  tken  we 


may  faind  n  complex,  number  xJP  x^,...p  x n  Mick  tkat  aQ  4  0  and 


1  n 

Pi(x)=axn  +  a.1xn~+...+a.  =  a  n  [x 
n  o  1  n  o  ■  , 

x=  1 


Xj  J 

A. 


If  P ( x )  satisfies  all  conditions  of  (2.7)  with  the  exception  of 
the  last,  i .e. , 

(m  -!) 

P  (xn)  =  0  , 


then  we  must  have 


m.  m,  mn  ,  m  -1 

P(x)  =  A(x)  (x  -  xQ)  (x  -  x-j )  1  ...  (x  -  xn-1)  (x  -  xn) 


where  A(x)  is  a  polynomial.  By  examining  the  degree  of  this  product, 
it  appears  that  A  is  a  constant.  Since,  moreover, 


(m  -1) 
p  n  (xn) 


A(mn  ‘  1)1  (xn  -  xo> 


m. 


(x. 


-  x  ,  ) 
n-1 ' 


m 


n-1 


=  0 


and  since  x.  f  x.,  i  f  j ,  we  have  A  =  0  and  therefore  P(x)  =  0.  The 

J 

homogeneous  interpolation  problem  has  the  zero  solution  only  and  hence 
the  nonhomogeneous  problem  possesses  a  unique  solution. 


2.3  Lagrangian  Interpolation  Polynomial 

This  polynomial  is  a  special  case  of  (2.2)  with  j  =  0.  Without 
any  ambiguity  the  superscript  o  on  <J>  and  the  second  subscript  o 
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on  y  in  (2.2)  may  be  deleted.  Thus  the  Lagrangian  interpolation 
problem  may  be  defined  by 


(2.8)  4>(xi)  =  y. 


for  i  =  0,1 ,. . . ,  n 


Consider  the  following  polynomial  of  degree  n, 


n 

n  (x  -  x.) 
i  =0  1 

(2.9) 

\(x) 

.  m. _ 

n 

n  (xk  -  x.) 
i=0  1 

W 

for  k  =  0,1 

n 


It  is  obvious  that 


(2.10) 


r 0 

if  k  i  j 

6kj 

if  k  =  j 

where  6^.  may  be  referred  to  as  the  Kronecker  delta  [11].  Therefore 
(2.8)  will  be  satisfied  if 

(2.11)  <D(x)  =  L  (x)  =  l  yk  £k(x) 

n  k=Q  K  K 


The  polynomial  L  (x)  may  be  called  the  Lagrangian  interpolation 
polynomial . 
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In  some  applications  of  interpolation,  the  tabular  points  may 
be  equally  spaced  and  in  this  case  (2.11)  can  be  simplified  [19]. 

Let  -  x.j_i  =  h  for  i  =  1,2 . n.  With  x  =  xq  +  th  (2.9)  may 

be  wri tten 


Ak(x) 


(-l)n_1  ,nN  t(t-l) _ (t-n) 

t-i  V  n! 


and  hence  (2.11)  becomes 


(2.12)  Ln(x)  -  (n*  l)(n+V  jo  (?)y. 


The  Lagrangian  formula  (2.11)  is  invariant  under  a  linear  trans¬ 
formation.  If  we  let  x  =  hu  +  a  and  x.  =  hu^  +  a  then  (2.9)  becomes 

n 


n  h(u  -  u. ) 
i=0  1 


(2.13) 


^ 


"  £k(u) 


n  h(uk  -  U.) 
i  =0  K  1 
i7k 


for  which  (2.11)  may  be  written  as 


Ln(x)  =  I  yk  Vx) 

n  k=0  K  K 


l  yk  Mu) 

k=0  K  K 


■  Ln<u> 


(2.14) 
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2.3.1  Uniqueness  of  the  Lagrangian  System  The  uniqueness  of 
the  solution  of  this  system  may  be  obtained  as  a  result  of  the  follow¬ 
ing  theorem,  the  proof  of  which  may  be  found  in  Davis  [11]: 

Theorem  2.4  Given  [n  +  1)  cLutlncl  points  x..,  I  =  0.1.....  n 

1  A* 

and  in  +  1 )  vainer  y.,  l  =  0,1 ,...,  n,  tken  them  exists  a  unique  poly- 

A* 

nomial  ( x )  oft  degree  n  fion  u)hlek 

(2.15)  P  [x . )  =  y.  I  =  0,1,...,  n 

Vi  A^ 

The  (n  +  1)  given  data  points  are  sufficient  to  specify  uniquely  a  poly 
nomial  of  degree  not  in  excess  of  n;  and  if  it  is  to  assume  the  values 

of  y.j  for  x  equal  to  x^ ,  i  =  0,1 . n,  its  coefficients  a^  for 

i  =  0,1,...,  n  must  be  such  that 


n  i 

(2.16)  l  a.  x'.  =  y.  for  j  =  0,1,...,  n 

i=0  1  J  J 

or  expressed  in  matrix  form, 


— 

— 

1 

x2 

-*S 

A  e  o 

0 

ao 

yo 

1  X] 

*1  •• 

ai 

(2.17) 

• 

= 

• 

1  xn 

x2 

xn 

a 

n 

n 

an 

yn 

-  .  i 

■ 
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where  the  determinant  of  the  coefficient  matrix  is  known  as  the 
Vandermonde  determinant. 

When  x.j ,  y.  for  i  =  0,1,...,  n  are  numeric,  (2.17)  becomes  a 
system  of  (n  +  1)  linear  equations  in  (n  +  1)  unknowns  and  can  be 
effectively  solved  by  Gaussian  elimination  with  scaling,  pivoting  and 
iterative  improvement  if  desired  [24].  If  all  the  x^'s  are  distinct, 
the  determinant  of  the  coefficient  matrix  in  (2.17)  is  nonsingular 
and  the  solution  a.  is  unique.  Since  the  interpolation  problem  (2.15) 
has  a  unique  solution,  no  matter  which  of  the  finite-difference  inter¬ 
polation  formulae,  e.g.  Newton's,  Gauss's,  Bessel's,  Stirling's,  etc., 
is  employed  for  a  given  set  of  points,  the  polynomial  used  will  be 
the  Lagrangian  polynomial. 

In  continuous  piecewise  interpolation  the  values  of  y.  .  in  (2.1) 

'  J 

are  not  given  for  j  >  1.  But  an  expression  may  be  obtained  for  these 
quantities  in  terms  of  y..  by  use  of  the  Lagrangian  interpolation 
formula  as  shown  in  the  next  chapter. 


CHAPTER  III 


ALGEBRAIC  DEVELOPMENT 


The  method  of  continuous  piecewise  interpolation  presented 
herein  is  an  extension  of  concepts  developed  by  Snyder  in  [32]  and 
[34].  The  basic  technique  is  continuous  transformation  of  a  segment 
of  a  polynomial  curve  of  an  arbitrary  chosen  degree,  as  interpolation 
proceeds  across  many  data  points  covering  the  full  range  of  data. 

The  transformation  of  the  segment  is  achieved  by  first  shifting  the 
segment  a  specified  number  of  points  to  the  right  and  then  re-evalua¬ 
ting  the  interpolation  formula.  The  derivatives  at  the  junction  point 
of  the  two  segments  are  equivalent  and  thus  the  process  is  continuous. 
In  addition,  the  application  of  symmetry  is  employed  in  the  derivation 
as  suggested  by  Charmonman  and  Wojtiw  [10]. 

3.1  Preliminary 

Through  each  of  the  subsets  of  (n  +  1)  points  of  the  set  of  data 
it  is  possible  to  pass  a  polynomial  of  degree  n.  The  number  of  points 
(n  +  1)  can  take  on  any  numeric  value.  For  convenience,  however,  the 
same  number  of  points  could  be  used  in  all  the  subsets  considered. 
Instead  of  requiring  the  Lagrangian  polynomial  Ln(x)  to  pass  through 
the  (n  +  1)  points  of  the  subset  we  may  use  the  Hermitian  polynomial 
Hn(x)  to  pass  through  (k  <  n)  points  for  k  >  2  and  to  satisfy  addi¬ 
tional  (n  +  1  -  k)  conditions  on  the  derivatives  at  the  first  and 
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last  of  the  k  points.  The  additional  (n  +  1  -  k)  conditions  on  the 
derivatives  are  not  given  explicitly  and  thus  must  be  obtained  by 
use  of  Lagrangian  polynomials.  Let  q  denote  the  number  of  points 
passed  through  by  the  Lagrangian  polynomial.  The  parameter  q  may 
be  either  odd  or  even,  but  for  convenience  will  be  taken  to  be  odd. 

Let  s  =  t  q/23  +  1.  Then  the  subset  of  k  equally  spaced  points  passed 
through  by  the  Hermitian  polynomial  may  be  defined 

(3.1)  (x. ,  i  =  s ,  (s  +  1 ),...,  (s  +  k  -  1 )  such  that  x.  >  x.  if  i  >  j} 

I  I  J 

such  that 

(3.2)  xs+i  =  x$+._1  +  h  i  =  1,2 .  (k  -  1) 


or 


(3.3) 


Vi  =  xs  +  1h  1  =  1)2 . (k  - 1J 


To  avoid  having  a  different  number  of  derivatives  specified  at  the 
first  point,  i.e.  x$ ,  than  at  the  last  point,  i.e.  x^^,  we  will 
require  that  if  (n  +  1)  is  odd  then  k  must  be  odd  and  if  (n  +  1)  is 
even  then  k  must  be  even.  Thus  in  either  case  (n  +  1  -  k)  will  be 
even.  For  simplicity  of  notation,  we  introduce 


m 


n  + 


1 

7 


(3.4) 


16 


and  define  the  order  of  the  multi  step  in  terms  of  the  subset 
of  k  points  i .e. , 

Definition  3.1  The.  ondeA  p  oft  muZtutzp  otf  any  fioAmula  may  be. 

denned  <x6 

(3.5)  p  =  k  -  1 

wkeJie.  k  tt>  the.  numbeA  ofi  point*  pa**ed  through  by  the.  HeAmltlan 
polynomial. 

Thus  we  use 


(3.6) 


q  =  n  +  1  -  p 


points  in  each  of  the  two  Lagrangian  polynomials  to  evaluate  the  m 
derivatives  specified  at  x$  and  at  x  -j.  Therefore  the  Hermitian 
polynomial  is  required  to  satisfy 


(3.7) 


^xs+i^  ~  ^s+i 


i  =  0,1,...,  (k  -  1) 


and 

(3.8)  4)^J^(xs+.)  =  y^j  for  i  =  0  and  (k  -  1) 

j  =  1 ,2, ... ,  m 
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One  Lagrangian  polynomial  is  first  passed  through  the  subset 


S]  =  (yi ,  i  =  0,1,...,  (q  -  1)} 


and  another  through  subset 


S2  =  i  =  (n  -  q  +  1),  (n  -  q  +  2),...,  (n  +  1)} 


while  the  Hermitian  polynomial  passes  through  ,  x  +p...* 


xs+k-l * 


3.1.1  Use  of  Symmetry  For  convenience  in  the  derivation  pro¬ 
cess  as  well  as  for  efficient  computation  of  the  coefficients  of 
the  polynomial  the  tabular  points  in  any  subset  of  (n  +  1)  points  may 
be  ordered  symmetrically.  If  the  number  of  points  is  even,  the 
subscripts  may  be 


(3.9) 


+  1  ,  .  .  .  , 


1 


2 


» •  •  • » 


and  if  the  number  of  points  is  odd,  the  subscripts  may  be 


(3.10) 


n 


"  7 


+  1 


» 


-1,  0,  1 


3.2  Laqrangian-Hermitian  Approach 

3.2.1  Sample  Derivation  The  use  of  Lagrangian  and  Hermitian 


. 
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interpolation  formulae  in  deriving  the  continuous  piecewise  inter¬ 
polation  formula  will  be  illustrated  by  an  example.  Consider  any 
set  of  six  points  symmetrically  ordered,  i.e.  y_2  5,  y_-|  5>  y_  5> 

y  5s  yl  5’  y2  5  tw0  subsets  si  =  (y -j »  i  =  -2.5  to  1.5)  and 

S2  =  (y. ,  i  =  -1.5  to  2.5)  as  in  Figure  3.1.  In  this  example  we 
require  continuity  of  the  first  and  second  derivatives  at  the  points 
x  =  -.5  and  x  =  .5.  Let  J-j(x)  represent  the  Lagrangian  polynomial 
passing  through  the  points  in  subset  S-j  and  J2(x)  the  Lagrangian 
polynomial  for  subset  S2. 


Figure  3.1  A  Set  of  Six  Points 
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Hence  the  two  Lagrangian  polynomials  of  degree  four  may  be  constructed 
by  use  of  (2.9)  and  (2.11)  with  =  i  for  i  =  -2.5  to  2.5. 


J 1  (x) 


(x  +  1 . 5) (x  +  . 5) (x  -  .5) (x  -  1.5) 

- r-i)c-2)T-3)(--4l - 


-2.5 


,  (x  +  2.5)  (x  +  .5)  (x  -  .5)  (x  -  1.5)  .. 

D(-i)  -2)(-3)  y 


-1.5 


(3.11) 


(x  +  2.5) (x  +  1 . 5) (x  -  ,5)(x  -  1.5) 

+  - (2)(l)(-l)(-2) -  y 


-.5 


(x  +  2 . 5) (x  +  1 . 5) (x  +  . 5) (x  -  1.5) 

+  - (3)(2)(1)(-1 ) - -  y 


.5 


and 


(x  +  2.5) (x  +  1 . 5 ) ( x  +  .5) (x  -  .5) 

+  - (4)(3)(2)(1) - -  y1.5 


I  f.,')  _  (x  +  . 5) (x  -  . 5) (x  -  1 . 5) (x  -  2.5)  .. 

J2(x)  "  - (-l)(-2)(-3)(-4)  y-l .! 


(x  +  1.5)(x  -  .5) (x  -  1 .5) (x  -  2.5) 

+  - (1)(-l)(-2)(-3) - y 


-.5 


(3.12) 


(x  +  1 .5) (x  +  . 5) (x  -  1 .5) (x  -  2.5) 

+  - (2}(l)(-1Jl-2) —  -  y 


.5 


(x  +  1 . 5 ) ( x  +  .5) (x  -  . 5) (x  -  2.5)  .. 

+  - (3)(2)(l)(-l) - -  y 


1.5 


(x  +  1 . 5) (x  +  .5) (x  -  . 5) (x  -  1.5)  .. 

+  - (4)(3)'(2llT) -  -  y2. 
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Using  (3.11)  we  obtain 

(3.13)  J-j  (-.5)  =  y 


(3.14)  Jj(-.5)  =  T7  (y-2.5  '  8y-1.5  +  8y.5  '  y1.5) 

and 


(3.15)  J-j  1  (- . 5)  -  T7  ("y-2.5 ,+  16y-l  .5  “  30y-.5  +  16y.5  "  yl  .5^* 

Similarly  (3.12)  gives 


(3.16) 

02( .5) 

(3.17) 

J£(.5) 

=  T?(y-1.5 

and 

(3.18) 

J£‘(.5) 

=  T?  (_y-l. 

=  y.5 

-  8y._5  +  8y15  -  y2  5) 

+  16y__5  -  30y  j  +  16y1>5  -  y2_5). 


The  second  set  of  equations  (3.16)  to  (3.18)  may  be  easily  obtained 

from  the  first  set  (3.13)  to  (3.15)  with  a  change  of 

indices  i.e.  y^  is  replaced  by  y^+-|  for  i  =  -2.5  through  to  1.5.  This 


. 

•'  s'  '  U<v  *  i*r  ’  £  vi  +,.|.t-)-r-  *  (•  .)'$&  (8f  E) 
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set  of  six  equations,  i.e.  (3.13)  through  to  (3.18),  may  be  represented 
in  matrix  form  as 


(3.19) 


1  - 

0,(-.5) 

0 

0 

12 

0 

0 

0 

y-2.5 

J{(-.5) 

1 

-8 

0 

8 

-1 

0 

y-l .  5 

J-j 1  (-.5) 

1 

-1 

16 

-30 

16 

-1 

0 

y-.5 

=  77 

J2(.5) 

0 

0 

0 

12 

0 

0 

y.5 

0£(.5) 

0 

1 

-  8 

0 

8 

-1 

y1.5 

jy(.5) 

0 

-1 

16 

-30 

16 

-1 

y2.5 

_ 



We  now  construct  a  Hermitian  polynomial  of  degree  five  which 
satisfies 


(3.20)  H^(-.5)  =  jjj)(-.5) 
and 

(3.21)  H*j)(.5)  =  Jp'h.5) 


where  j  =  0,1,2  . 
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One  way  to  construct  H^(x)  is  to  consider  a  Lagrangian  interpolation 
polynomial  L^(x)  which  assumes  values  y  y  3  at  the  points 
x_  5,  x  g.  Therefore  applying  (2.9)  and  (2.11)  L-j(x)  may  be 
written  as 

(3.22)  L,(x)  =  ^(-2x  +  1)  y__5  +  ^(2x  +  1)  y_5  . 

Obviously  the  difference  H^(x)  -  L-j(x)  is  a  polynomial  of  degree  five 
which  will  vanish  at  the  points  x_  5  and  x  Thus 

(3.23)  H5(x)  -  L.|  (x)  =  (x  -  x_  5)(x  -  x  5)  H3(x) 

where  H3(x)  is  a  Hermitian  polynomial  of  degree  three  to  be  deter¬ 
mined.  Hence 

(3.24)  H5(x)  =  ^(-2x  +  1)  y_>5  +  ^(2x  +  1)  y_5 

+  (x  +  .5)(x  -  .5)  H3(x)  . 

Differentiating  (3.24)  with  respect  to  x  we  have 

(3.25)  H£(x)  =  -y_  5  +  y  5  +  (x  +  .5)(x  -  .5)  H^(x) 


+  2x  H3(x) 


The  condition  in  (3.20)  with  j  =  1  is  obtained  by  equating  (3.25) 
wi th  x  =  - . 5  to  (3.14)  i . e . 
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(3.26) 


+  y.5  "  M’*5)  • 


Similarly  using  (3.21)  with  j  =  1,  equation  (3.25)  with  x  =  .5,  and 
(3.17)  yields 


(3.27) 


+  y.5  +  H3^ ,5) 


Hence  we  may  write 

(3.28)  H3(-.5)  =  y^“y-2.5  +  8y-1.5  "  12y-.5  +  4y.5  +  y1.5^ 

and 

(3*29)  H3(*5)  =  lVy-1.5  +  4y-.5  "  12y.5  +  8yl . 5  _  y2.5^  * 

Differentiating  (3.25)  with  respect  to  x  again,  substituting  x  =  -.5 
and  x  =  .5,  and  using  (3.20)  and  (3.21)  respectively  we  now  get 
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(3.30) 

H3^",5}  ~  ^“y-2.5  +  6y-.5  ~  8y.5  +  3y  1 . 5 ^ 

and 


(3.31) 

H^(.5)  -  ^-3y.1.5  +  8y._5-6y-5  +  y2_5)  . 

Thus  we  have  found  H^( x )  satisfying  (3.28)  through  to  (3.31).  If 
we  use  (2.9)  and  (2.11)  again  we  have 


(3.32) 

h3(x)  =  £<-2x  +  1)  H3(-.5)  +  j(2x  +  1)  H3(.5) 

+  (x  +  .5)(x  -  .5)  H-j  (x) 

where  H^(x)  is  a  Hermitian  polynomial  of  degree  one  to  be  determined. 
After  differentiating  with  respect  to  x  and  substituting  x  =  -.5  and 
x  =  .5  we  can  obtain  from  (3.32) 


(3.33) 

V"*5)  =  ^-(3y.2.5  "  14y-1.5  +  26y-.5  "  24y.5 

+  ]1yi.5  -  2y2.5} 

and 


(3.34) 

H,(.5)  =  7T<-2y.2.5  +  11y.l.5  •  24y-.5  +  26y.5 

-  14^1 .5  +  3*2.5>  • 
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Hence  H^x)  =  H]  (-.5)  +  (.5)  -  H]  (-.5))(x  +  .5) 

and  by  back-substitution  we  obtain  H^( x ) ,  which  enables  us  to 
obtain  ( x )  i .e. 

h5(x)  =  7^[x5(-160y_2  5  +  800y_1  5  -  1600y^  g  +  1600y  g 

-  80°y1>5  +  16°y2.5)  +  x406y_2.5  ~  48y-i.5  +  32y-.5  +  32y.s 

-  48y1>5  +  16y2  5)  +  x3(144y_2  g  -  848y_1 +  1824y_>5  -  1824y>5 

(3.35)  +  848-j  5  -  144y2  5)  +  x2(-40y_2>5  +  312y_K5  -  272y_  5  -  272y>5 

+  ^  1 2y  15  -  40y2  +  x(-26y_2  ^  +  162y_-|  ^  -  1124y_  ^ 

+  1124y.5  -  162*1.5  +  26y2.5]  +  (9y-2.5  -  75*-1.5  +  450*-5 
+  450y_5  -  75,. g  +  9y2  5)]. 

With  a.  defined  as  the  coefficient  of  the  Hermitian  polynomial  of  degree 
five 

5 

H5(x)  =  l  a.  x1 

0  i=0  1 


(3.36) 


(3.35)  may  be  expressed  in  matrix  form 


9 

-  75 

450 

450 

-  75 

-  26 

162 

-1124 

1124 

-162 

(3-37)  75S 

-  40 

144 

312 

-  272 

-  272 

312 

-848 

1824 

-1824 

848 

16 

-  48 

32 

32 

-  48 

. 

-160 

800 

-1600 

1600 

-800 

3.2.2  General 

Deri 

vation 

We  now 

cons 

as  stated  by 

(3.7) 

and  (3.8)  subject  to 

the 

specified  in 

section  3.1 

.  The 

Lagrangi 

an  i 

9 

26 
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' 

“““ 

y-2.5 

ao 

y-l  .5 

al 

y-.5 

a2 

y.5 

a3 

y1.5 

a4 

y2.5 

a5 

may  be  expressed  in  symmetric  form.  For  the  number  of  points  (n  +  1) 
we  may  write 


(3.38) 


n/2 

•n(x)  =  l  y,*  Ux) 

n  i=-n/2  1  1 


where 


. 
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Also  if  J-|(x)  denotes  the  Lagrangian  polynomial  -|(x)  passing  through 
the  points  in  subset  S-j  =  {y^. ,  for  the  first  q  points  of  (n  +  1)}  and 
J2 (x )  for  the  subset  =  (y^ *  for  the  last  q  points  of  (n  +  1)}  then 
for  the  subset  of  k  points  we  may  evaluate  k  equations  i.e. 


(3.39)  J-j(x)|  >  J  (x)|  J  (x )  I  x-x  »  *MX^  I  x~x 

!  X  Xs  W  X-Xs+1  W  X-X$+k_2  2  X-X$+k_i 


where  w  may  be  either  1  or  2  depending  on  which  is  more  appropriate. 
This  corresponds  to  the  left-hand  side  of  (3.7).  We  next  evaluate  the 
m  derivatives  of  (3.38)  at  the  two  points  xg  and  x  +k_-j,  by  first 
differentiating  the  equation  with  respect  to  x  and  then  evaluating 
the  required  derivative  at  these  two  points  to  yield  m  equations  for 
the  point  x$ ,  i.e. 


(3.40) 


j;(x)|x=x  ,  J]'(x)|x=x  ,...,  J™(x)|x=x 
s  s  s 


and  an  additional  m  equations  for  i.e, 


(3.41) 


Jo (x) |  ,  Jo'(x) L_y  Jo(x)| 

1  x_xs+k-l  L  x_xs+k-l  6  xs+k-l 


A  total  of  2m  equations,  i.e.  (3.40)  and  (3.41),  correspond  to  the 
left-hand  side  of  (3.8).  Thus  a  total  of  (2m  +  k)  equations  are  formed 
by  (3.39),  (3.40)  and  (3.41). 


•  !  . 


■ 

.  " 
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We  now  use  Lagrangian  formula  (3.38),  to  define  another  polynomial 
L^_ i ( x )  where  k  is  the  number  of  points  in  the  subset  of  k  points. 

We  further  assume  that  there  exists  a  polynomial  of  degree  n,  denoted 
by  H  (x)  satisfying  the  right-hand  side  conditions  of  equations  (3.7) 
and  (3.8).  The  difference  Hn (x )  -  -|(x)  will  be  a  polynomial  of 

degree  n,  which  will  vanish  at  the  points  x$+^ ,  i  =  0,1 .  (k  -  1). 

Consequently, 


Hn(x)  -  *- k_  1  ( x )  ■  wk(x)  Hn_k(x) 


where 


(3.42) 


For  any  polynomial  H  . (x)  the  function 


Hn(x>  =  Lk-l(x)  +  wk(x)  Hn-k(x) 


(3.43) 


will  assume  the  values  of  y^  at  the  points  of  interpolation.  We 

now  obtain  H  . (x)  such  that  the  remaining  conditions  are  satisfied, 
n-K 

by  first  differentiating  both  sides  of  (3.43)  with  respect  to  x  and 


evaluating  x  =  x^  where  i  =  s  and  (s  +  k  -  1)  i.e. 


(3.44) 
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but  k(x.)  =  0  since  =  0*  Since  )  /  0  we  find 

H  i,(x.)  at  each  point  where  H'(x.)  is  defined.  This  process  may 

be  continued  m  times  for  the  m  derivatives  of  the  right-hand  side  of 

(3.8).  Each  time  the  coefficient  of  the  highest  derivative  of 

Hn_k(x)  at  the  points  will  be  oo^x^).  Thus  the  problem  of  finding 

H  (x)  becomes  that  of  finding  H  . (x)  satisfying  the  conditions 

n  *  k 


(3.45) 


,  .  x  d j  ,Hn(x)  ~  Lk-l(x\ 

Hn-k  ^xs+i  ^  =  dxJ  Vx^ 


j  =  1,2,...,  (m  -  1 ) 

i  =  0  and  (k  -  1 ) 


and 


(3.46) 


Hn-k<xs+i>  ■ 


H  (x)  -  L,  ,(x) 
nv  '  k-lv  ' 


wk(x) 


i  =  0,1,...,  ( k  -  1 ) 


Similarly  the  system  (3.45)  and  (3.46)  may  be  reduced.  After  carry¬ 
ing  out  this  process  m  times,  it  is  possible  by  back-substitution  to 
obtain  the  final  result  in  a  matrix  form. 

3.3  Direct  Derivation  by  Use  of  Polynomial  Expressions 

A  more  direct  approach  to  obtain  a  continuous  piecewise  inter- 

1L 

polation  formula  is  to  assume  a  (v  =  q-lr n  degree  polynomial  Qy(x) 
passing  through  the  q  points  in  either  subset  S-|  or  i.e. 


Q  (x)  =  l  a.  x1 

v  i=0  1 


(3.47) 
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such  that  x  =  -g,  (-g  +  1 ) ,. . . ,  0,1 ,  g  where  g  =  j  and  q  has  been 
taken  for  convenience  to  be  odd.  Thus  for  subset 


(3.48)  q  (_g  +  i )  =  y  i  =  0,1 . v 

V  -  \  +1+P 

while  for  subset  S-j  (3.48)  holds  with  p  =  0.  Equation  (3.47)  may  be 
written  in  matrix  form  as 


(3.49) 


Qv(-g) 

Qv(-g+D 

• 

• 

Qv(°) 

• 

= 

0 

Qv(g-D 

JO 

< 

1 

(-g+D1 


0 


(g-n1 

(g)1 


(-g)2....(-g)v  ’ 

(-g+l)2..(-g+l)v 


0  . ..  0 

(g-i )2° . . (g-i  )v 
(g)2  »..(g)v 


a 


1 
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xs  for  subset  S-j  corresponds  to  x  =  0  in  (3.47).  We  may  also  calculate 
the  inverse  of  the  coefficient  matrix  in  (3.49)  to  obtain  the  coef¬ 
ficient  aQ  in  (3.47)  which  represents  the  left-hand  side  of  (3.7)  in  terms 
of  y -'s  at  x  .  Coefficients  a.  after  multiplication  by  the  factorial 

*  <j 

of  the  i ^  power  of  x  in  (3.47)  yield  the  left-hand  side  of  (3.8)  in 
terms  of  y. ' s  at  xc .  Similarly  we  can  let  x  =  0  correspond  to  the 
xs+k-l  Point  in  the  subset  and  obtain  the  same  coefficients  as  for 

x  ,  except  for  y  n  now  replaced  by  y  n  .  .  Following  this  procedure 

"7  ‘7  P 

for  (3.7)  and  (3.8)  we  may  obtain  the  system 


(3.50)  c  =  Dy 

where  c  is  a  vector  of  (2m  +  k)  conditions  as  given  by  the  left-hand 
side  of  (3.7)  and  (3.8),  D  a  matrix  of  dimension  (n  +  l)-by-(2m  +  k) 
expressed  in  terms  of  the  derived  coefficients,  and  y  a  vector  of 
length  (n  +  1)  of  the  y. 's.  This  is  equivalent  to  passing  a  Lagrangian 
polynomial  through  the  points  in  subsets  S-j  and  S^. 

We  now  assume  that  another  polynomial, 

2m+k-l 

(3-51)  H2m+k-l(x)  '  J0  ei  * 

provides  the  basic  form  for  the  Hermitian  interpolation  function.  The 
(2m  +  k)  coefficients  of  this  polynomial  satisfy  the  requirements  given 
by  (3.7)  and  (3.8).  We  obtain  k  equations  for  (3.7)  and  2m  equations 
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for  (3.8).  The  total  of  (2m  +  k)  equations  may  be  written  as 

(3.52)  d  =  Bz 

where  d  is  a  vector  of  the  (2m  +  k)  conditions  as  given  by  (3.7)  and 
(3.8),  B  a  matrix  of  order  (2m  +  k)  expressed  in  terms  of  the  derived 
coefficients,  and  z  =  {eQ,  e-j,...,  e2m+k- 1  * 

Since  c  =  d,  from  (3.7)  and  (3.8),  then  (3.50)  and  (3.52)  with 
F  =  B’^D  becomes 

(3.53)  z  =  Fy 

where  z  =  {eo»  e^ . e2m+k-l^T*  F  a  matnx  dimension  (n  +  l)-by- 

(2m  +  k)  giving  numeric  coefficients,  and  y  the  vector  of  y. 's  of 
length  (n  +  1).  Hence  (3.53)  is  the  matrix  equation  defining  the  various 
interpolation  formulae  and  F  is  its  coefficient  matrix.  Thus  (3.51)  can 
now  be  expressed  in  terms  of  the  tabular  points  y-j's.  It  can  be  noted 
that  the  coefficient  matrix  F  will  be  square  if 

(3.54)  q  =  2m+l  , 
and  that  a  unique  solution  exists  only  if 


(3.55) 


2m+k  <  n+1 . 
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3.3.1  Sample  Derivation  If  we  now  consider  the  same  example 
as  in  sub-section  3.2.1,  then  (3.49)  may  be  written 


(3.56) 


y-2.5 

1 

-2 

4 

-8 

16  ~ 

ao 

y-1.5 

1 

-1 

1 

-1 

1 

al 

y- ,  5 

1 

0 

0 

0 

0 

a2 

y .  5 

1 

1 

1 

1 

1 

a3 

y1.5 

1 

2 

4 

8 

16 

J 

a4 

— 

for  the  subset  S-j  =  (y. ,  i  =  -2.5  to  1.5}.  The  expression  for 
S>2  =  {y^ ,  i  =  -1.5  to  2.5}  is  exactly  the  same  as  (3.56)  except  for 
the  left-hand  side  where  y^  is  replaced  by  y^  for  i  =  -2.5  to  1.5. 
The  inverse  of  (3.56)  becomes 


(3.57) 


Thus 


ao 

0 

0 

24 

0 

0 

y-2.5 

al 

1 

=  2T 

2 

-16 

0 

16 

-2 

y-1.5 

a2  j 

-1 

16 

-30 

16 

-1 

y-.5 

a3 

-2 

4 

0 

-  4 

2 

y .  5 

a4 

1 

4 

6 

-  4 

1 

y1.5 

[_ 

mm, 

_ 

<(,(-.5)  =  (°y-2.5  +  0y-1.5  +  24y-.5  +  0y.5  +  0y1.5^ 
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while 

<t>'(-.5)  =  ^  (2y_2  5  -  16y,1>5  +  0y_>5  +  16y_5  -  2y,  g)  *  1! 

12  (y-2.5  “  8y-l . 5  +  8y .  5  '  y1.5^ 

and 

♦  "f--5)  =  2T  ('y-2.5  +  16y-1.5  "  30y-.5  +  16y.5  ’  y1.5*  *  2! 

=  T7  (-y-2.5  +  1 6y- 1 .5  ‘  30y- . 5  +  16y.5  ‘  y1.5> 

Similarly,  4>(.5),  4> '  ( . 5 )  and  <J>''(.5)  may  be  obtained  to  yield  the 
coefficient  matrix 


-1  16  -30  16  -  1 

1  -  8  0  8  -1 


0 

0 


(3.58) 


D  = 


T? 


0 

0 


0  12  0  0  0 
0  0  12  0  0 


0  1_8  0 
0  -  1  16  -30 


8  -1 

16  -1 


For  the  basic  form  of  the  interpolation  function  we  assume  a  general 
fifth-degree  polynomial 


H5(x) 


X 


i 


(3.59) 


/ 
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for  which  we  may  obtain  the  matrix  form 


(3.60) 


—  - 

— 

r* 

H5"(-.5) 

0 

0 

64 

-96 

96 

-80 

H5'(-.5) 

0 

32 

-32 

24 

-16 

10 

H5(-.5) 

1 

=  12 

32 

-16 

8 

-  4 

2 

-  1 

H5(.5) 

32 

16 

8 

4 

2 

1 

H5'(.5) 

0 

32 

32 

24 

16 

10 

H5"  (.5) 

0 

0 

64 

96 

96 

80 

- 

e2 

e3 

e4 

e5 


The  coefficient-matrix  of  (3.60)  is  B  in  (3.52).  We  can  now  calculate 
the  inverse  of  B  and, with  D  from  (3.58),  obtain  equation  (3.37). 


3.4  The  Computer  Algorithm 

A  computer  program  may  be  written  using  the  approach  presented 
in  section  3.3  to  obtain  the  coefficient  matrix  F  of  (3.53)  in  floating¬ 
point  representation. 

Algorithm  3.1  Given  the  three  variables  k,  q,  and  m  where 
k  is  the  number  of  points  passed  through  by  the  Hermitian  polynomial, 
q  the  number  of  points  used  in  evaluating  the  derivatives,  and  m 
the  highest  order  of  the  derivative  required  for  continuity. 

Step  1:  Test  if  the  selection  of  k,  q,  and  m  yields  a  unique  solution. 
If  not,  then  exit. 

Step  2:  Calculate  the  coefficient  matrix  B  in  (3.52),  which  provides 
the  basic  form  for  the  Hermitian  interpolation  function,  as 
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Set  up  matrix  B  of  order  (2m  +  k)  containing  zeros 
and  calculate  x  =  i  -  (k  +  1),  i  =  1,2,...,  k 
i.e.  the  tabular  points  through  which  the  Hermitian 
polynomial  passes. 

Evaluate  x^  in  (3.51)  for  j  =  0,1,...,  (2m  +  k  -  1) 

i.  L 

using  x  in  step  2.1  and  place  these  in  the  (m  +  i)tn 
row  of  B.  This  gives  the  requirements  of  equation 
(3.7). 

Evaluate  the  m  derivatives  at  x  =  -  A(  k  -  1 )  and 
x  =  Jjr  (k  -  1)  and  place  these  in  the  first  m  and 
the  last  m  rows  of  B  respectively.  Thus  the 
requirements  of  (3.8)  are  satisfied. 

Step  3:  Now  calculate  the  coefficient  matrix  D  in  (3.50).  This  may  be 
considered  as  passing  the  Lagrangian  polynomial  through 
subsets  S-j  and  S ^ 

Step  3.1:  Set  up  the  coefficient  matrix  D  of  dimension 
(k  +  q  -  l)-by-(2m  +  k)  such  that  D  =  { d .  . } , 

'  J 

j  =  1,2,...,  (k  +  2m)  and  i  =  1,2,...,  (k  +  q  -  1) 
yields  d. .  =  1  if  j  =  i  -  [m  +  i(l  -  q)]  and 

'  J  ^ 

d..  =  0  otherwise.  This  is  equivalent  to  satisfy- 
i  J 

ing  conditions  (3.7). 

Step  3.2:  Calculate  x  =  i  -  ^  (q  +  1 )  for  i  =  1,2,...,  q 

and  evaluate  xJ  in  (3.47)  for  j  =  0,1,...,  (q  -  1). 


follows: 
Step  2.1: 


S tsp  2.2: 


Step  2.3: 
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This  now  can  be  expressed  as  the  coefficient  matrix 
of  (3.49). 

Step  3.3:  Obtain  the  inverse  of  this  coefficient  matrix  by 
Gaussian  elimination. 

J.  L. 

Step  3.4:  Use  the  (i  +  1 ) Lri  rows  of  this  inverse  times  i! 

for  i  =  1,2,...,  m  to  obtain  the  m  derivatives  of 
(3.8)  at  x  =  -  -^(k  -  1)  and  at  x  =  |-(k  -  1),  and  place 
these  in  the  first  m  and  the  last  m  rows  of  D 
respectively.  Hence  the  coefficient  matrix  D  is 
set  up. 

Step  4:  Calculate  B~^  by  Gaussian  elimination  and  obtain  F  =  B’^  D. 

F  is  now  the  desired  coefficient  matrix  in  floating-point  form. 

A  flow  diagram  for  such  a  program  is  given  in  Figure  3.2. 

3.5  Discussion  of  Error  in  Machine  Derivation 

We  now  investigate  the  computational  errors  arising  in  the 
derivation  of  a  continuous  piecewise  interpolation  formula  by  a 
digital  computer.  All  the  computations  and  results  presented  herein 
were  performed  on  an  IBM  S/360/67  computer  using  APIA 360  or  Iverson's 
language  [14],  [20].  In  APL\360  all  floating-point  computations  are 
performed  in  extended  long  word  (64  bits  with  56  bits  of  fractional 
part)  and  the  final  result  may  be  rounded  and  displayed  in  1  to  17 
significant  figures.  The  results  presented  herein  on  error  bounds 


t 
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Figure  3.2  Flow  diagram  for  Algorithm  3.1  for  the 

computation  of  the  coefficient  matrix  F  of 
(3.53)  in  floating-point. 
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are  due  to  Wilkinson  [36]  and  Forsythe  and  Moler  [15]. 


3.5.1  Calculation  of  Matrix  D  Since  g  and  i  in  (3.47)  are 
integers  the  coefficient  matrix  W  of  (3.49)  (step  3.2) 


W  = 


(-g) 

(-g+i) 


(g-D 

(g) 


(-g) 

(-g+i )' 


(g-i V 
(g)2 


(-g)v 
(-g+i )v 


(g-i  )v 
(g)v 


may  be  obtained  without  any  error  if  a  special  routine  using  integer 
arithmetic  is  used.  We  now  compute  the  inverse  of  W  by  using  Gaussian 
elimination  [24].  In  this  method  the  scaled  matrix  W  is  decomposed 
into  the  product  of  two  triangular  matrices  L  and  U.  Pivoting  only 
involves  permutations  of  the  row  subscripts.  Since  permutation  is 
irrelevant  to  the  error  analysis  we  may  simply  neglect  pivoting  in 
our  analysis.  The  decomposition  consists  of  computing  a  sequence 
of  matrices  W^  =  W,  W^,  ...,  W^,  where  the  matrix  W^  has 
zero  elements  below  the  diagonal  in  the  first  (k  -  1)  columns.  The 
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(k+1)  (k) 

matrix  Wv  '  is  obtained  from  Wv  by  subtracting  a  multiple  of  the 
k^*1  row  from  each  of  the  rows  below  it  and  the  rest  of  is  left 
unchanged.  The  multipliers  are  chosen  so  that  if  there  were  no  round¬ 
ing  errors  would  have  zeros  below  the  diagonal  in  the  k**1  column. 

Let  =  {wjkb  and 


(3.62) 


mi.k  =  fl  (w1\k/wM)  > 


i  >  k  +  1 


r 


0 


for  i  >  k  +  1 ,  j  =  k 


(3.63)  w?kt^  =  (  fl  (w!k|  -  m.  .x  wbb  for  i  >  k  +  1  and  j  >  k  +  1 

1»J  1»J  K,J 


(k) 

Vwij 


otherwise 


Also  let 

(3.64)  U  =  W(n) 


and 
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1 

(3.65) 

L  = 

m2,1 

m3,l 

1 

m3,2 

1 

m  , 
n,l 

m  o  . . 

. . . .  1 

n,2 

* 

Obviously  L  and  U  are  lower  and  upper  triangular,  respectively.  The 
following  theorem  is  stated  and  proven  by  Forsythe  and  Moler  [15]. 

Theorem  3. 1  The  Pollution  x  computed  by  Gaussian  elimination 
with.  pivoting  saXijs{ie s  tkc  equation 

(3.66)  (W  +  6W)  x  =  b  . 

furthermore 

(3.67)  ||  6W  ||  <  I.OHn3  +  3 n2)  p  ||W  ||  „  u  . 

where 

p  *  max  |wjfeJ.|/  lid'll  „ 
i.j.k  *->1 

and  u  the  unit  round- i.e. 

1  j_  + 

u  =  y  3  (bounded  operations ) 

u  =  3* ~ ^  [chopped  operations  ) 
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On  the  IBM  S/360/67  using  APL\36Q,  the  value  of  3  -  2  and  t  =  53.  For 
even  approximate  equality  to  hold  in  (3.67)  the  "worst"  round-off  must 
occur  at  each  step  of  the  calculation  and,  furthermore,  equations 
like  ||  AB||  ^  =  ||  A||  •  ||  B||  ^  must  hold.  Forsythe  and  Moler  state 

that  no  example  is  known  where  |[  6 W ||  ^  is  even  close  to  the  given 
bound.  Wilkinson  [36,  p.  108]  states  that  ||6W||  ^  is  rarely  larger 
than  nu  ||  W||  oo.  He  further  gives  the  bound  of  6W  for  sufficiently 
large  n  as 

(3.68)  ||  6W||  ^  <  gi  2  (2.005n2  +  n3  +  ,25n4  2  *’ )  j|  U||  ^ 

,  (kl  i 

where  g-j  is  the  largest  element  of  any  |WV  '  |  and 

"tl  t 

(3.69)  2  1  =  (1.06)  2’1 

Also  he  remarks  that  the  last  term  in  the  parenthesis  of  (3.68)  is 

-t  ”^1  3 

unimportant  if  n2  «  1,  and  the  dominant  term  is  g-j  2  n  .  The 

dominant  term  comes  from  errors  in  the  solution  of  the  triangular 

system  of  equations.  The  bound  in  (3.68)  may  only  be  attained  in 

exceptional  circumstances,  so  that  even  the  replacement  of  the  factor 
2  3 

2.005n  +  n  by  its  square  root  gives  a  somewhat  higher  value  than 
would  normally  be  expected  in  practice.  Either  the  results  of  Wilkinson 
or  Forsythe  and  Moler  may  be  used  to  give  a  bound  for  the  calculation 
of  W"1.  Let  E-j  be  the  error  in  obtaining  W"  .  The  coefficient  matrix 
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D  obtained  by  step  3.4  in  Algorithm  3.1  of  section  3.4  will  not  result 
in  any  additional  errors.  Thus  the  error  in  obtaining  the  matrix  D 
is  also  E-j . 

3.5.2  Calculation  of  Matrix  F  =  ET^  D  The  matrix  B  in  step  2 
of  Algorithm  3.1  of  section  3.4  can  also  be  obtained  without  any  error. 

But  its  inverse  computed  by  Gaussian  elimination  may  contain  an  error, 
denoted  by  Eg,  for  which  the  theorem  in  sub-section  3.5.1  will  also 
hold.  A  final  error,  denoted  by  Eg,  may  arise  when  the  products  of 
two  matrices  B“^  D  is  calculated.  Thus  we  have 


(3.70) 


fKB-1  D)  e  B"1  D  +  E 

=  (B_1  +  E2)(D  +  E-j )  +  E3 

=  B  1  D  +  EgD  +  B  ^E-j  +  Eg 


Thus  the  total  error  E  is 


(3.71) 


E  =  EgD  +  B  ^  E-j  +  Eg 


Let  us  consider  E-|  =  Eg  =  0.  This  is  equivalent  to  assuming  that 
no  error  arises  in  calculating  the  inverse,  since  it  is  possible  to 
obtain  the  inverse  of  an  integer  matrix  in  integer  form.  This  is  done 
in  hand  computation  by  performing  all  calculations  to  the  common 
denominator.  Then  (3.70)  becomes 
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(3.72)  fHB"1  D)  =  B~!  D  +  £3 

where  E3  is  the  error  due  to  the  product  of  B~^  D.  For  the  product 
of  two  matrices  in  floating-point  Wilkinson  [36]  gives  the  result 

(3.73)  ||  E3||  E  <  2  n  ||  B"1 1|  E ||  D 1 1  £ 

where  t-j  satisfies  (3.69)  and  the  Euclidean  norm  for  any  matrix  A  is 

(3.74)  II  A||  c  =  (  l  l  |a..|2)1/2 

1  j  ^ 

The  bounds  given  are  hardly  likely  to  be  attained  since  ||  ||  £||  D||  £ 

may  be  very  much  greater  than  1 1  B“ ^  D 1 1  E  and  the  bound  for  E3  may  be 
much  greater  than  the  actual  error. 


3.6  Accuracy  of  Machine  Derivation 

Although  no  reasonable  bound  of  round-off  error  can  be  obtained 
we  can  check  the  actual  error  of  the  result  of  the  machine  derivation. 
It  is  possible  to  follow  Algorithm  3.1  given  in  section  3.4  using 
hand  computation.  Such  computations  may  be  carried  out  in  integer 
arithmetic  and  expressed  in  terms  of  a  common  denominator.  Thus 
the  results  obtained  are  without  any  round-off  error  at  all.  Compari¬ 
son  of  the  exact  result  to  the  floating-point  machine  derivation 
revealed  that  the  actual  errors  in  the  floating-point  result  were  so 
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small  that  it  is  possible  to  convert  the  floating  point  result  into 
integer  form  by  multiplying  it  by  a  certain  number,  a. 

If  a  is  the  multiplying  factor  such  that  the  floating-point  repre¬ 
sentation  may  be  replaced  by  an  integer  form  then  it  is  possible  to 
obtain  a  numerical  value  for 

(3.75)  max  |aa..  -  b.  .  |  <  e..  all  i,  j 

•  ♦  I  .1  I  .1  “  I  .1 


and 

(3.76)  $  =  min  e.. 

1 3 

where  F  =  {a.  .}  and  G  =  (b.  .}  i.e.  F  is  the  computed  coefficient  matrix 
*  J  '  J 

in  floating-point  while  G  is  the  exact  one  obtained  by  hand  computation 
in  integer  form.  The  results  of  the  cases  considered  for  3  and 
||  aF  -  G 1 1  ^  (see  (3.74))  are  given  in  Table  III. 2  of  Appendix  III. 
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CHAPTER  IV 


CONTINUOUS  PIECEWISE  INTERPOLATION  FORMULAE 

Various  continuous  piecewise  interpolation  formulae  can  be  con¬ 
structed  either  in  integer  or  floating-point  notation.  The  formulae 
depend  upon  three  basic  parameters,  k,  q,  and  m,  where  k  is  the  number 
of  points  through  which  the  Hermitian  polynomial  passes,  q  the  number 
of  points  used  in  evaluating  the  derivatives,  and  m  the  highest  order 
of  the  derivative  required  for  continuity.  Other  quantities,  i.e. 
total  number  of  points,  the  order  of  the  multistep  etc.,  may  be  used 
but  these  may  be  expressed  in  terms  of  k,  q,  and  m.  Thus  the  various 
formulae  presented  may  be  classified  according  to  the  format 

(4.1)  (k,  q,  m}  . 

In  order  that  a  unique  solution  will  exist  the  parameters  are  subject 
to  the  conditions  given  in  Chapter  III.  The  formulae  presented  in 
this  chapter  are  for  k  up  to  five,  q  up  to  eleven  and  m  up  to  three. 
Other  formulae  with  higher  values  of  the  three  parameters  may  be  ob¬ 
tained  in  a  similar  manner. 

For  identification  purpose  any  formula  of  the  mLn  kind  is  defined 

th 

to  be  the  formula  which  has  continuity  up  to  the  rn  n  derivative.  The 
formulae  of  the  m  kind  may  be  sub-divided  according  to  the 
order  of  the  multi  step.  For  classifying  the  order  (k  -  1)  of  the 
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multistep, adjectives  such  as  single,  double,  triple,  etc.  will  be  used 
to  modify  the  noun  "formula".  Further  modification  of  the  form 
(n  +  1  =  k  +  q  -  l)-po1nt  will  be  used  to  indicate  the  total  number  of 
points  used  in  the  subset.  As  an  example,  a  formula  requiring  continuity 
up  to  the  second  derivative,  with  the  order  of  the  multistep  being  one 
and  using  six  points  would  be  called  the  six-point  single-step  formula 
of  the  second  kind.  The  results  for  the  cases  considered  are  summarized 
in  Table  III.l  of  Appendix  III. 

The  use  of  symmetry  as  mentioned  in  Subsection  3.1.1  results  in 
the  absolute  value  of  the  coefficient  matrix  being  symmetric  in  a  sense 
not  the  same  as  in  matrix  algebra  (see  Appendix  II).  If  the  total 
number  of  points  used  is  even,  symmetry  is  about  a  line  between  the 
centre  most  columns.  For  an  odd  number  of  points  the  symmetry  is 
about  the  middle  column. 

Also  for  each  of  the  continuous  piecewise  interpolation  formulae 
it  is  possible  to  generate  a  set  of  desk-calculator  tables  on  an  elec¬ 
tronic  computer.  A  procedure  for  obtaining  such  tables  is  given  in 
section  4.4  of  this  chapter.  An  example  of  such  tables  may  be  found 
in  Table  III. 3  of  Appendix  III. 

4.1  Continuous  First  Derivative 

A  number  of  different  continuous  piecewise  interpolation  formulae 
with  continuous  first  derivative  may  be  obtained.  As  mentioned  before, 
these  may  be  called  formulae  of  the  first  kind  since  m  =  1.  The 
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simplest  possible  formula  for  odd  q  is  when  k  =  2  and  q  =  3.  Hence 
a  total  ofn+l=k+q-l=4  points  (x_-j  g,  x_  g,  x  g,  x-j  g}  are 
used  with  the  first  derivative  at  x  g  expressed  in  terms  of  the  first 
three  points  and  at  x  g  in  terms  of  the  last  three  points.  The 
Hermitian  polynomial  satisfies  four  conditions  and  passes  through  the 
two  interior  points.  By  format  (4.1)  this  formula  may  be  called 
Formula  2.3.1  or  the  four-point  single-step  formula  of  the  first  kind. 
Here  the  interpolating  polynomial  is  moved  to  the  right  one  data  point 
at  a  time  accounting  for  the  phrase  "single-step".  The  coefficient 
matrix  F  of  (3.53)  in  this  case  is 


-1 

9 

9 

-1 

1 

2 

-22 

22 

-2 

16 

4 

-  4 

-  4 

4 

-8 

24 

-24 

8 

and  is  square  as  may  be  expected  from  (3.54).  This  formula  is  exact 
to  a  polynomial  of  degree  two  since  the  lowest  degree  of  the  polynomials 
used  is  two. 

As  an  extension  of  Formula  2.3.1  we  may  consider  q  =  5  for  which 
Formula  2.5.1  or  the  six-point  single-step  formula  of  the  first  kind, 
may  be  obtained.  As  in  the  case  of  Formula  2.3.1  the  interpolating 
polynomial  also  moves  one  data  point  at  a  time  to  the  right.  Formula 
2.5.1  is  exact  to  a  polynomial  of  degree  three  since  the  lowest  degree 
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of  the  polynomials  used  is  three.  In  this  case  the  coefficient  matrix 
F  of  (3.53)  is  rectangular, 


1 

-  9 

56 

56 

-  9 

1 

1 

-2 

14 

-128 

128 

-14 

2 

96 

36 

-  32 

-  32 

36 

-4 

8 

-56 

128 

-128 

56 

-8 

Similarly  we  may  obtain  coefficient  matrices  for  q  =  7,  9  and 
11  and  denote  these  as  Formulae  2.7.1,  2.9.1,  and  2.11.1  respectively. 

All  these  formulae  are  single-step  and  exact  only  to  a  polynomial  of 
degree  three. 

Another  extension  of  Formula  2.3.1  is  to  increase  the  order  of 
multistep(p  =  k  -  1).  Let  us  consider  p  =  2,  i.e.  k  =  3,  then  the 
interpolating  polynomial  is  now  moved  two  tabular  points  at  a  time 
to  the  right  instead  of  one  as  for  p  =  1.  The  simplest  case  possible 
for  odd  q  is  q  =  5  corresponding  to  Formula  3.5.1  or  the  seven-point 
double-step  formula  of  the  first  kind.  Here  the  Hermitian  polynomial 
passes  through  three  points  instead  of  two  as  in  the  single-step  formula. 
For  symmetry  the  seven  points  used  may  be  ordered  as  (x^,  x_^,  x_p 
X  x-j ,  X2»  x3).  The  formula  thus  obtained  is  exact  to  a  polynomial 
of  degree  four  with  a  coefficient  matrix  F  in  (3.53)  being 
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0 

0 

0 

48 

0 

0 

0 

-i 

8 

-37 

0 

37 

-8 

1 

(4.4) 

1 

i 

-8 

47 

-80 

47 

-8 

1 

i 

-8 

13 

0 

-13 

8 

-1 

-i 

8 

-23 

32 

-23 

8 

-1 

For  all  double-step  formulae  the  total  number  of  points  will  be  odd. 
Similarly  as  for  single-step  formula  we  may  increase  the  value  of 
q  to  obtain  formulae  which  will  exact  to  a  polynomial  of  degree  four 
but  which  use  more  points  in  evaluating  the  derivative.  Likewise  it 
is  possible  to  extend  the  double-step  formulae  such  that  the  inter¬ 
polating  polynomial  moves  three  tabular  points  at  a  time  to  the  right, 
and  these  formulae  may  be  called  the  triple-step  formulae  of  the  first 
kind.  With  q  =  5  we  have  the  eight-point  triple-step  formula  of  the 
first  kind,  or  Formula  4.5.1,  with  the  coefficient  matrix 

— 

-  27  216  -  918  4185  4185  -  918  216  -  27 

18  -144  756  -8586  8586  -  756  144  -  18 

120  -960  3888  -3048  -3048  3888  -960  120 

-  80  640  -3232  7056  -7056  3232  -640  80 

-  48  384  -  864  528  528  -  864  384  -  48 

32  -256  832  -1440  1440  -  832  256  -  32 


<4- 5) 
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The  eight  points  may  be  ordered  {x^ ,  i  =  -3.5  to  3,5}  where  subset 
(x.. ,  i  =  -3.5  to  .5}  is  used  in  evaluating  the  derivative  at  x  -j  ^ 
while  the  subset  {x. ,  i  =  -.5  to  3.5}  is  used  in  evaluating  the  de¬ 
rivative  at  x-j  The  Hermitian  polynomial  passes  through 
x. ,  i  =  -1.59  -.5,  .5,  1.5.  This  formula  is  exact  to  a  polynomial 
of  degree  four. 

Two  examples  of  quadruple-step  formulae  of  the  first  kind. 
Formulae  5.3.1  and  5.5.1  may  be  found  in  Table  II. 1  of  Appendix  II. 

4.2  Continuous  First  and  Second  Derivatives 

Formulae  with  continuity  of  the  first  two  derivatives  may  be 
called  the  formulae  of  the  second  kind.  With  k  =  2,  the  simplest 
case  possible  with  odd  q  is  Formula  2.5.2  which  was  expressed  in  non- 
symmetric  form  by  Snyder  [34]  and  subsequently  in  symmetric  form  by 
Charmonman  and  Wojtiw  [10].  This  formula  may  be  called  the  six-point 
single-step  formula  of  the  second  kind  with  coefficient  matrix 


9 

-  75 

450 

450 

-  75 

9 

-  26 

162 

-1124 

1124 

-162 

26 

1 

-  40 

312 

-  272 

-  272 

312 

-  40 

76F 

144 

-848 

1824 

-1824 

848 

-144 

16 

-  48 

32 

32 

-  48 

16 

-160 

800 

-1600 

1600 

-800 

160 

. 
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This  formula  is  exact  to  a  polynomial  of  degree  four. 

Like  the  formulae  of  the  first  kind,  the  formulae  of  the  second 
kind  may  be  extended  by  increasing  values  of  q  and  k.  An  example  of 
such  an  extension  is  the  eight-point  triple-step  formulae  of  the  second 
kind  or  Formula  4.5.2,  which  has  the  coefficient  matrix 


(4.7) 


1 

21)736 


-162 

1053 

-  3483 

12960 

12960 

-  3483 

1053 

-162 

135 

-  918 

3564 

-27783 

27783 

-  3564 

918 

-135 

756 

-4860 

15228 

-11124 

-11124 

15228 

-4860 

756 

-636 

4296 

-16032 

31068 

-31068 

16032 

-4296 

636 

-448 

2672 

-  5328 

3104 

3104 

-  5328 

2672 

-448 

400 

-2592 

7360 

-11920 

11920 

-  7360 

2592 

-400 

64 

-  320 

576 

-  320 

-  320 

576 

-  320 

64 

-  64 

384 

-  1024 

1600 

-  1600 

1024 

-  384 

64 

The  interpolating  polynomial  here  moves  three  tabular  points  at  a  time, 
while  the  Hermitian  polynomial  passes  through  x^ ,  i  =  -1.5,  -.5,  .5,  1.5. 
This  formula  is  exact  to  a  polynomial  of  degree  four. 

4.3  Continuous  First  Three  Derivatives 

Formulae  with  continuity  up  to  the  third  derivative  may  be  called 
formulae  of  the  third  kind.  Three  examples  of  such  formulae  are 
given  in  Table  II. 3  of  Appendix  II,  one  single-step  formula  and  two 
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double-step  formulae.  Formula  2.7.3,  is  single-step  and  exact  to  a 
polynomial  of  degree  six,  while  Formulae  3.7.3  and  3.9.3,  are  double¬ 
step  and  exact  to  a  polynomial  of  degree  six  and  eight  respectively. 

4.4  Generation  of  Desk-Calculator  Tables 

Tables  of  continuous  piecewise  interpolation  coefficients  for 
desk-calculator  use  may  be  obtained  on  an  electronic  computer.  Consider 
the  set 

R  =  {Xj ,  j  =  -  y*  ‘  7  +  1 » •  •  • »  such  that  xj  =  1  i f  j  =  -  £  +  1  and 
x.  =  0  otherwise}  where  i  is  one  of  the  values  0,1,...,  n.  For  the 

J 

coefficient  matrix  F  in  (3.53)  the  individual  and  unit  contribution 

of  y  n  .  for  an  appropriate  x  is  obtained  by  considering  the  set  R 
-7+ 1 

and  solving  equation  (3.53)  for  all  values  of  e,  i.e. 

J 


2m+k-l 

C.(x)  =  l  e.  xJ 

1  j=0  J 


(4.8) 


The  quantity  c..(x)  is  called  the  interpolating  coefficient  [34].  Since 
there  are  (n  +  1)  of  these  unit  contributions  we  may  interpolate  in 


any  set  {x.. ,  i  =  -  j. 


. . . , 


by 


n 


(4.9) 


■mo 

■ 

q**9 
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Desk-calculator  tables  can  be  generated  in  any  increment  of  x 
from  x  =  -.5  to  x  =  .5  for  single-step  formulae;  from  x  =  -1.0  to 
x  -  1.0  for  double-step  formulae;  from  x  =  -1 .5  to  x  =  1 .5  for  triple¬ 
step  formulae  and  from  x  =  -2.0  to  x  =  2.0  for  the  quadruple-step 
formulae.  An  example  of  coefficients,  for  Formula  2.5.2,  as  generated 
on  a  computer  is  given  in  Table  III. 3  of  Appendix  III. 


CHAPTER  V 


TEST  RESULTS  AND  CONCLUSIONS 

A  practical  test  of  the  method  of  continuous  piecewise  interpola- 
tion  formulae  is  difficult  to  design.  In  practical  applications  the 
true  function  is  not  known,  since  if  it  were  there  would  be  little 
point  in  having  an  interpolation  formula  because  any  desired  values 
could  be  computed  from  the  true  function.  Thus  a  test  to  decide  which 
continuous  piecewise  interpolation  formula  would  give  the  best  approxi¬ 
mation  to  intermediate  values  for  measurements  of  physical  quantities 
would  be  impossible  to  construct.  Therefore  the  various  continuous 
piecewise  interpolation  formulae  are  tested  on  known  functions. 

5 . 1  Test  Functions  and  Procedure 

The  known  functions  considered  were 

(5.1)  y  =  x1  for  1  =  6,7,10 


and 

(5.2)  y  =  ex 

For  (5.1)  interpolation  was  carried  out  in  the  interval  [1.025,2]  in 
increments  of  .025,  while  for  (5.2)  in  the  intervals  [10.025,11]  and 
[25.025,  26]  also  in  increments  of  .025. 
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The  number  of  multiplications  required  in  the  application  of  the 
formula  can  be  minimized  by  using  the  mirror- image  property  of  the  coef¬ 
ficient  matrix  F  in  (3.53).  In  other  words,  with  F  =  {f.  .}  denoting  the 

■  J 

coefficient  matrix, 

n+1 

(5.3)  e.  =  l  f..  y  1-1.2 . (2m+k) 

1  j-1  1J  -£+j-l 


and  symmetry  in  F,  (5.3)  may  be  written 


(5.4) 


t-1 


f  l  ^(y  n  ±  yn  ) 
j-l  10  4+  j-l  7+l'j 


ei 


f? 

t-1 


fityo+  l  fii<y  n  1  yn  > 

0  j-l  1J  -f  +  j-l  £  +  1-j 


for  (n+1)  even 


for  (n+1 )  odd 


n+1 

where  t  =  [-£— ]J+  1  with  plus  and  minus  signs  depending  upon  the  coef¬ 
ficient  matrix  F.  The  term  f..  y  in  (5.4)  is  zero  for  i  even  and  (n+1) 

I  V  U 

odd.  Also  the  expression 


t-1 


I  f-,(y  n  ±  y  ) 

j=i  1J  4  +  j-i  ^  +  i-j 


for  (n+1)  odd  in  (5.4)  is  zero  for  i  =  1.  Thus  for  (n+1)  even  (5.4)  re¬ 
quires  (n+l)+2  multiplications  and  n  additions  (algebraic)  as  opposed  to 
(n+1)  multiplications  and  n  additions  in  (5.3).  For  (n+1)  odd,  in  the 
worst  case  (5.4)  requires  (-y4+l  multiplications  and  n  additions  as 
opposed  to  (n+1)  multiplications  and  n  additions  in  (5.3).  Therefore 
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computation  of  the  e.'s  by  use  of  (5.4)  can  be  made  to  be  about  twice 
as  efficient  as  in  using  (5.3).  Furthermore  a  much  more  efficient  way 
to  evaluate  any  polynomial  is  to  factor  out  the  x's.  For  example,  for 


(5.5) 


Y(x)  =  l  a.x 
isl 


i-1 


we  may  write 


(5.6)  Y(x)  =  a-j  +  x(a2  +  x(a^  +  x(a^  +  x(a&+  xa^)) )) . 

(5.6)  requires  only  five  multiplications  and  five  additions  as  opposed 
5 

to  l  i  =  15  multiplications  and  five  additions  for  (5.5).  In  some 
i=l 

machines  where  one  multiplication  takes  much  longer  than  one  addition 
then  addition  time  can  be  neglected  and  the  computation  of  Y(x)  in  (5.6) 
requires  about  1/3  of  that  for  (5.5).  In  some  other  machines  one  multi¬ 
plication  takes  about  as  long  as  one  addition  and  the  computation  of 
Y(x)  in  (5.6)  costs  about  half  that  in  computing  Y(x)  in  (5.5).  The 
number  of  multiplications  required  for  computing  each  formula  by  the 
above  procedures  are  given  in  Table  IV. 1  of  Appendix  IV. 

The  test  procedure  was  to  apply  the  various  formulae  to  the  given 
function  and  obtain  interpolated  values  for  this  function  in  increments 
of  .025.  Since  we  know  the  true  solution,  we  may  then  obtain  the  rela¬ 
tive  error  for  each  interpolated  value.  We  now  calculate  the  largest 
relative  error  for  each  formula  and  use  this  as  a  measure  of  how  good 
an  approximation  has  been  obtained.  The  results  of  (5.1)  are  given  in 
Table  IV. 2  of  Appendix  IV,  while  those  of  (5.2)  in  Table  IV. 3  of  the 
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same  appendix. 

5.2  Conclusions 

The  following  observations  may  be  drawn  from  the  results  of  this 
study: 

1.  Comparison  of  the  exact  result  to  the  floating-point  machine 
derivation  revealed  that  the  actual  errors  in  the  floating-point  result 
were  so  small  that  it  is  possible  to  convert  the  floating-point  result 
into  exact  integer  form. 

2.  Many  new  formulae  presented  in  this  thesis  gave  a  better  approxi¬ 
mation  than  Snyder's  Formula  2.3.1  [32],  but  these  required  a  larger 
number  of  multiplications.  For  example.  Formula  2.5.1  gave  a  more  accu¬ 
rate  result  for  the  test  examples  than  Formula  2.3.1  in  both  Tables  IV. 2 
and  IV. 3  of  Appendix  IV.  However  Formula  2.5.1  requires  15  multiplication 
as  opposed  to  11  for  Formula  2.3.1.  Similarly  for  Snyder's  Formula  2.5.2 
[34],  there  are  cases  where  a  more  accurate  approximation  is  achieved  re¬ 
quiring  a  larger  number  of  multiplications  i.e.  Formula  2.7.2  in  Tables 
IV. 2  and  IV. 3.  Formula  3.7.1  requires  the  same  number  of  multiplications 
and  yet  it  gives  a  more  accurate  result  than  Formula  2.5.2  for  the  test 
examples.  No  formula  was  found  requiring  a  lower  number  of  points  and 
giving  a  better  approximation  than  Formula  2.5.2. 

3.  Increasing  the  number  of  points  in  evaluating  the  derivative 
i.e.  q,  produces  a  better  approximation  up  to  a  certain  value  of  q,  after 
which  no  improvement  results  in  the  approximation.  This  is  obvious  if 

we  consider  Formulae  2.3.1,  2.5.1,  2.7.1,  2.9.1,  and  2.11.1  for  x^  and  yJ 
in  Table  IV. 2,  and  ex  in  Table  IV. 3.  The  same  degree  of  accuracy  may 
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be  achieved  using  either  Formula  2.5.1  or  Formulae  2.7.1 s  2.9.1,  or  2.11.1. 
Thus  when  more  than  one  formula  give  the  same  accuracy  we  should  use  the 
formula  requiring  the  smallest  number  of  multiplications. 

4.  Formulae  of  the  second  kind  gave  a  lower  maximum  relative  error 
than  those  of  the  first  kind,  as  is  obvious  for  Formulae  2.7.1  and  2.7.2 
in  Tables  IV. 2  and  IV. 3.  Comparing  formulae  of  the  second  kind  with  those 
of  the  third  kind,  no  noticeable  improvement  resulted  in  the  maximum  re¬ 
lative  error,  but  an  increase  in  the  number  of  multiplications  was  observed. 
This  would  imply  that  formulae  of  the  second  kind  can  be  used  instead  of 
the  first  kind  to  obtain  a  better  approximation,  and  instead  of  the  third 
kind  to  obtain  a  decrease  in  the  number  of  multiplications. 

5.  An  increase  in  the  order  of  the  multistep  for  constant  values  of 

q  and  m  does  not  necessarily  produce  an  improvement  in  the  maximum  relative 
error.  An  improvement  may  result  depending  upon  whether  p  is  even  or  odd, 
the  value  of  q,  and  the  value  of  i  in  x1 .  If  we  consider  the  maximum 
relative  errors  for  Formulae  2.3.1,  3.3.1,  4.3.1,  and  5.3.1,  for  an  in¬ 
crease  in  p  from  one  to  three  no  improvement  results,  while  for  an  in¬ 
crease  in  p  from  3  to  4  an  improvement  is  observed  for  Table  IV. 2  only. 

If  on  the  other  hand  we  consider  Formulae  2.5.1  and  3.5.1  an  improvement 
results  in  the  maximum  relative  error  with  an  increase  in  the  order  of 
the  multistep. 

6.  We  may  obtain  the  formulae  requiring  the  lowest  number  of  multi¬ 
plications  for  the  lowest  maximum  relative  error  for  Tables  IV. 2  and  IV. 3 
of  Appendix  IV.  These  for  Table  IV. 2  for  x6,  x7,  and  x10  are  Formulae 
3.7.2,  3.9.3,  and  3.9.2  respectively,  while  for  Table  IV. 3  in  the  intervals 
[10.025,11]  and  [25.025,26]  the  formulae  are  Formulae  3.9.3  and  2.7.2. 
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5.3  Suggestions  for  Future  Research 

1.  A  computer  program, say  in  Fortran  IV,  may  be  written  such  that 
all  operations  are  performed  in  integer  notation.  This  would  mean  that 
all  operations  would  have  to  be  carried  out  under  the  lowest  common 
denominator.  This  would  result  in  the  coefficient  matrix  being  in  integer 
form.  Special  care  must  be  given  to  the  problem  of  overflow  in  many 
parts  of  the  program. 

2.  The  various  continuous  piecewise  interpolation  formulae  may  be 
tested  on  measurements  of  physical  guantities  such  as  streamflow,  pre¬ 
cipitation,  etc.  For  any  particular  engineering  problem  we  may  make 
the  measurements  at  some  intermediate  point  in  addition  to  the  egui- 
distance  set  of  points.  Then  test  the  various  formulae  available  to  find 
the  formula  which  reguires  the  smallest  number  of  multiplication  from 
among  those  giving  results  with  an  acceptable  degree  of  accuracy. 
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APPENDIX  I 


APLX360  FUNCTIONS 

The  following  pages  contain  a  description  and  a  listing  of  APL\360 
functions  used  throughout  this  thesis.  Two  functions  SIGDIG  and  EFT  are 
omitted  since  they  are  found  in  LIBRARY  1.  The  function  SIGDIG  controls 
the  number  of  significant  digits,  between  one  and  seventeen,  in  the  final 
result,  while  the  function  EFT  enables  us  to  output  the  results  with  a 
different  number  of  significant  digits  in  each  column.  For  example,  if 
we  have  four  columns  we  may  output  the  first  column  say  with  two  significant 
digits  while  the  other  three  with  eight  significant  digits  as  in  Table 
III .3  of  Appendix  III . 

CPI  B  corresponds  to  Algorithm  3.1  of  subsection  3.4  in  Chapter  III 
where  B  is  the  vector  of  the  three  parameters  in  (4.1).  This  function  uses 
INVERSE  M  to  calculate  the  inverse  of  a  non-singular  square  matrix  M  by 
Gaussion  elimination. 

A  VALABETA  B  and  A  EANORMACPIF  B  give  the  values  of  beta  in  (3.83) 
and  ||  otF-G 1 1 ^  respectively ,  A  is  the  vector  of  the  three  parameters  in  (4.1) 
while  B  is  the  reciprocal  of  the  coefficient  of  the  formula.  This  coef¬ 
ficient  is  given  for  each  formula  in  Appendix  II. 

The  function  C  DESKACOEFF  B  is  used  in  generating  the  desk-coefficients 
of  Formula  2.5.2  given  in  Table  III. 3  of  Appendix  III.  C  is  a  vector  of 
three  components  representing  the  starting  point,  the  number  of  rows,  and 
the  group  of  three  columns  desired  in  the  table  of  coefficients.  B  is  the 
vector  of  the  three  parameters  in  (4.1).  This  function  calls  Z  POLY  T  to 
calculate  the  value  of  a  polynomial  with  its  coefficients  represented  by 


' 
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the  vector  T  at  the  point  Z. 

The  function  BAPPLICACPIF  A  is  used  in  obtaining  the  values  of 
Tables  IV. 2  and  IV. 3  (i.e.  the  results  of  (5.1)  and  (5.2)  respectively). 

B  is  a  vector  of  three  components  representing  the  first  point  of  inter- 

1.  L. 

polations  the  increment,  and  the  l  power  of  x  for  Table  IV. 2.  For  the 
results  of  Table  IV. 3  the  third  component  of  B  is  replaced  by  ©  x  (for 
exp  x)  in  the  function  and  deleted  from  B  in  the  header.  A  is  the  vector 
of  the  three  parameters  in  (4.1).  The  largest  relative  error  for  each 
formula  is  obtained  in  the  output. 
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V  APPLIC  ACPI FLU ] V 


V 


[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
[9] 
CIO] 
[11] 
[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 

[32] 


Q+B  APPLICLCPIF  A \C \D  ;E  ;F  ;G \H  ;I  ;«7 \K \L \M ;N \0  ;P \R \S \T  ;U 

iViW;XiYiZ 

J+-\-  ((P/1)t3) 

C+-1  5  p  0 

Q+(J,  1  )pP+K+M«-N*-0 
IJ+CPI  V<~A  [  P  +  i  3  ] 

^CB[1]+B[2]xii0x(4-(7C1]=4)))*B[3] 
f/+(-(0.5  +  ((7[l]-2)x(  (P«-0 .25)+0.25))))+ 

0. lxiiox(7[i]-i) 

Z^(S[l]-((L(7[2]v2))+l)x5[2]xl0)+(3[2]xi0xi((7[l]-l)+ 

VZ21)) 

D+-Z  *B  [  3  ] 

E+2  |((p£/)[2]) 

F+  L  ((pi/)[l])-s-2 

G«-L  (  (p£/)[2]  )*2 

->  (  E=  0  )  / 1  5 

ClHl+UZH+i  (2xl(F+l))-l);i(G+l)]+.x(((Z?[iG]+Z?[l  +  <|>(i<7)  + 
Gl)  )  ,0[G+1]  ) 

•>17  9CiE]+UiH+2x  \Fi  iG]+.  x(Z?[  i  G  ]  -Di  1  +  4>  (  x  G  )  +  <7  ]  ) 

C[tf  >£/[#«-(  (2xiF)-l  )  ;iG]  +  .x(£>[  iG]+P[<|>(  iG)+G]) 
CLHl+UZH+2x\Fi  iG]  +  .  x(Z?[  i  G] -£>[<!>(  iG)+G]  ) 

0+-J7  POLY  C 

->(  (Y«-L4*(y[l]-1  )  )  =  4  )  /  2  0 
->2  1  ,I«-1 
T*-l+F-*-0 
->(  J  =  J)/30 

L+(Z+ 0  .  2  5  +F  +  Z  )  *5  [  3  ] 

->(P=0)/26 

C[P  >£/[#«-(  (2xx(F+1))-1);i(G+1)]  +  .x(((L[iG]+L[1  +  4>(iG)  + 
G]  )  )  , L [ G+l ]  ) 

->28  ,C[/?>P[P«-2x  iF;  XG]+  .  x  (L[  iG]-L[  1  +  <K  iG)+G] ) 

G[  #>£/[#«-(  (2xxF)-i);iG']+.x(LClfl!]+L[<J)(\G)+G]) 
G[ff>tf[ff«-2xiF;  iG]  +  .  x(L[  XG] -£[«>(  iG)+G]  ) 

POLY  C 

■>21,  J*-J  + 1 

£[M+-M+1  ;  1  XFC27]  XCT+Fi  T  /F«- 1  S+-(  *-G  )  ] 

P«-P+3 

+{J>K+K+1 ) /4 


[1] 


VPOLY [□ ] V 
V  G+-Z  P(7LY  2* 

X/(  ( (pZ)  .pDp^xZo 


*  1  +  x  pT 
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VCTimiV 

V  Q+CPI  BiK;M;N;DiC;A;V 

[1]  +( ( (2xB[3] )+S[l] )<( (B[l]+B[2] )-l ) )/3 

[2]  +OtpQ+('NO  UNIQUE  SOLUTION  EXISTS ') 

[3]  A+(KtK+(N+BZ  1]  )  +  2xD*-S[3]  )pO 

[4]  AID+  iN  ‘,l  +  ((\N)-0  .  5xN+l  )  °  .*V+~1+\K 

C5]  Al  \  D;  >4>[1]  ("l*(  -iD)°  .  +V)*A[D  +  N+\D  ;  X  (§(KtD  )p  !  \D)x(  \ 
D)°.lV)x(0.5*N-l)*(-\D)o.+V 

[6]  C+(iK)o .=( iN+M-1 )+D+0 . 5x1 -M+B[ 2] 

[7]  C[  iD  ;iM]+<p[l]C[D+N+iD  iN-l-iM]<-((S}(M,D)p  !  x D  )  x ( INVERSE ( ( 
\M)-0. 5xM+l  )o  .  *”l+  \M)  [  1+  iZ? ;  ] 

[8]  Q+-  ( INVERSE  A)  +  .xC 

V 


V  IN  VERSE l D]V 

V  R+INVERSE  M\N\A\K\I\V 

[1]  Y«-(  ,MZ  iA+Nl+A  o  .=/!),(  .M4-MC  ji4,i4])  ,(  pAf)  C*>1  ]  ) 

[2]  -►(  (1£“30)<  |A/C^;Z]  +  Oxp(  (  ,#[*■, I ;  >M[  J,£  ;  ]  )  ,  (  tI+(K-l)  +  V\ 
t/V+\ML (K-l )  +  \N+l-K;ia  ) ) )/4 

[3]  +0 ,pR+( ' DETERMINANT  IS  ZERO'  ) 

[4]  Y«-Ox((  ;KlxK*A)o  .xMUC-,1)  ,(  tMLK  ;]+MtK  iK]  ) 

) 

[5]  +(NZK+K+1  )/2 

[6]  7?+A/[  ;4+tf] 

V 


V  7i4L  A5£T/1  [  □  ]  V 

V  VALABETA  B  \  U  \  V 

[1]  Q+UIV\\ /V+\U+%  ( ( L  0 . 5  +  {CPI  A) xB ) - ( ( CPI  A)xB))l 

V 


VELN0RMLCPIFZU1V 

V  Q+A  E LNORMACPIF  B\C 

[1]  C<r%  ( (  L  0 . 5+  (CPI  A)  xB  )  -  (  (  CPI  A)xB)) 
C2]  Q+(+/(CxC))* 0.5 

V 


'  •  ■ 

. 

'  • .  f  f 
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VDESKACOFFFZUlV 

V  Q*-C  DESKACOEFF  B  ;A  ;  Z  ;  CO  ;  Cl  ;  C2 

[1]  A^CPI  B 

[2]  Z«-( -CZ  1 ] )  +  0 . Olx \CZ 2] 

[3]  C 0«-Z  POLY  A [  ;C[ 3 ]  +  1  ] 

[4]  Cl+Z  POLY  AZ  iCZ 3]+2] 

[5]  C2+-Z  POLY  AZ  ;CC3]  +  3] 

[6]  ('  *  C'iCZ3l+li' 

3 ]  +  2  ;  '  C 1  ;C[3]  +  3) 

[7]  e+(10  2  ,6p  16  8  )EFT§(  4 , £[ 2 ]  )  p  Z  ,  CO  t  Cl  t  02 

V 


C'  \Cl 


, 
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APPENDIX  II 


COEFFICIENT  MATRICES  OF  THE  FORMULAE 

The  various  continuous  piecewise  interpolation  formulae  are  classi¬ 
fied  according  to  the  highest  order  of  the  derivative  required  for  con- 
tinuity  and  are  called  formulae  of  the  mL  kind.  The  coefficient  matrix 
F  of  (3.53)  for  each  formula  is  given  in  integer  form  in  this  appendix. 
The  formulae  of  the  first  kind  are  given  in  Table  II. 1,  the  formulae  of 
the  second  kind  in  Table  II. 2,  and  the  formulae  of  the  third  kind  in 
Table  II. 3. 
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TABLE  II .1 

FORMULAE  OF  THE  1st  KIND 


FORMULA  2,3.1 

COEFFICIENT  OF  MATRIX  1*16 


"l 

9 

9 

"1 

2  “22 

22 

"2 

4 

4 

"4 

4 

"8 

24 

"24 

8 

FORMULA 

2.5.1 

COEFFICIENT 

OF 

MATRIX 

1 

*96 

1 

"9 

56 

56 

"9 

1 

"2 

14 

128 

128 

"14 

2 

"4 

36 

"32 

"32 

36 

"4 

8 

"56 

128 

"128 

56 

"8 

FORMULA  2.7.1 

COEFFICIENT  OF  MATRIX  1*480 


1 

10 

54 

285 

285 

54 

10 

1 

2 

16 

72 

630 

630 

72 

16 

"2 

4 

"40 

216 

180 

180 

216 

40 

4 

8 

64 

288 

600 

"600 

288 

"64 

8 

TABLE  II.l  (continued) 
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TABLE  II. 1  (continued) 

FORMULA  3.3.1 

COEFFICIENT  OF  MATRIX  1*8 

0  0  8  0  0 

1  “6  0  6  "1 

"1  8  ~14  8  "1 

~1  2  0  ”2  1 

1  “4  6  "4  1 


FORMULA  3.5.1 

COEFFICIENT  OF  MATRIX  1*48 


0 

0 

0 

48 

0 

0 

0 

“1 

8 

~  3  7 
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37 

"8 

1 

1 

“8 

47 

"80 

47 

"8 

1 

1 

~  8 

13 

0 

13 
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"1 

"1 
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32 

"2  3 
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FORMULA 

3  .  7 

.  1 

COEFFICIENT 

OF 

MA  TRIX 

1 

*240 

0  0 

0 

0 

240 

0 

0 

0 
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1  "9 

46 

189 

0 

189 
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"1 
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CD 
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44 

231 

390 

231 

44 

9 

"1 

1  9 

46 

69 

0 

"69 

46 

"9 

1 

1  "9 

44 

111 

150 

"111 

44 

"9 

1 
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TABLE  II. 1  (continued) 


FORMULA  4  •  J  •  1 

COEFFICIENT  OF  MATRIX  1*3456 


81 

"459 

2106 

2106 

"459 

81 

"54 

378 

“4320 

4320 

"378 

54 

360 

1944 

"1584 

"1584 

1944 

"360 

240 

1616 

3648 

"3648 

1616 

"240 

144 

"432 

288 

288 

"432 

144 

"96 

416 

"768 

768 

"416 

96 

FORMULA  ^  .  d  .  i 

COEFFICIENT  OF  MATRIX  1*6912 


"27 

216 

"918 

4185 

4185 

"918 

216 

"27 

18 

144 

756 

8586 

8586 

"756 

144 

"18 

120 

960 

3888 

"3048 

3048 

3888 

"960 

120 

80 

640 

"3232 

7056 

7056 

32  32 

"640 

80 

48 

384 

"864 

528 

528 

"864 

384 

"48 

32 

"256 

832 

1440 

1440 

"832 

256 

"32 

FORMULA  °  °  •  1 

COEFFICIENT  OF  MATRIX  1*576 


0 

0 

0 

576 

0 

0 

0 

24 

136 

"488 

0 

488 

136 

24 

12 

"80 

500 

864 

500 

"80 

12 

30 

"158 

226 

0 

226 

158 

"30 

15 

94 

"241 

324 

"241 

94 

"15 

"6 

22 

26 

0 

26 

22 

6 

3 

"14 

29 

36 

29 

"14 

3 
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TABLE  II. 1  (continued) 


FORMULA 


5.5.1 


COEFFICIENT  OF  MATRIX  1*1152 


0 

0 

0 

0 

1152 

0 

0 

0 

0 

8 

”64 

272 

"960 

0 

960 

”272 

64 

”8 

"4 

32 

”160 

992 

”1720 

992 

”160 

32 

”4 

10 

80 

"316 

432 

0 

”432 

316 

”80 

10 

5 

”40 

188 

”472 

638 

”472 

188 

”40 

5 

2 

”16 

44 

”48 

0 

48 

”44 

16 

”2 

"1 

8 

”28 

56 

”70 

56 

”28 

8 

”1 
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TABLE  1 1. 2 

FORMULAE  OF  THE  2nd  KIND 


FORMULA 

2.5.2 

COEFFICIENT 

OF  MATRIX 

1*768 

9 

75 

450 

450 

'75 

9 

“26 

162 

"1124 

1124 

"162 

26 

"40 

312 

272 

272 

312 

"40 

144 

848 

1824 

"1824 

848 

"144 

16 

48 

32 

32 

"48 

16 

160 

800 

"1600 

1600 

"800 

160 

FORMULA 

2.7.2 

COEFFICIENT 

OF  MATRIX 

1*11520 

28 

275 

"1377 

6890 

6890 

"1377 

275 

"28 

80 

"614 

2430 

"16300 

16300 

2430 

614 

"80 

128 

"1240 

5832 

4720 

"4720 

5832 

"1240 

128 

"448 

3376 

1  252  8 

23840 

"23840 

12  528 

"3376 

448 

64 

560 

"1296 

800 

800 

1296 

560 

"64 

512 

3680 

11232 

"18880 

18  880 

11232 

3680 

"512 

FORMULA 

3  . 

5 . 2 

COEFFICIENT 

OF 

MA TRIX 

1 

t96 

0 

0 

0 

96 

0 

0 

0 

"3 

20 

"79 

0 

79 

20 

3 

4 

28 

124 

"200 

124 

28 

4 

4 

24 

36 
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"36 

24 

"4 

"6 

40 

106 

144 

106 

40 

6 

1 

4 

5 
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5 

4 

1 

2 

12 

30 

40 

30 

"12 

2 
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TABLE  II. 2  (continued) 
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TABLE  1 1. 3 

FORMULAE  OF  THE  3rd  KIND 


FORMULA  2.7.3 

COEFFICIENT  OF  MATRIX  1*92160 


225 

2205 

11025 

55125 

55125 

11025 

2205 

225 

818 

"6158 

23178 

136630 

136630 

"23178 

6158 

"818 

1036 

"9980 

46764 

"37820 

37820 

46764 

"99  80 

1036 

"5720 

41960 

145080 

265480 

"265480 

145080 

"41960 

5720 

"560 

4720 

"10800 

6640 

6640 

"10800 

4720 

"560 

12640 

"89248 

269280 

"450080 

450080 

"269280 

89248 

"12640 

64 

320 

576 

320 

320 

576 

"320 

64 

11392 

79744 

239232 

398720 

"398720 

239232 

"79744 

11392 

FORMULA  3.7.3 

COEFFICIENT  OF  MATRIX  1*11520 


0 

0 

0 

0 

11520 

0 

0 

0 

0 

83 

"690 

2890 

9802 

0 

9802 

"2890 

690 

"83 

139 

1176 

475  6 

16424 

25410 

16424 

4756 

1176 

"139 

125 

990 

"3670 

4870 

0 

4870 

3670 

990 

125 

269 

"2232 

8492 

"18184 

23310 

18184 

8492 

"2232 

269 

49 

342 

878 

926 

0 

926 

"878 

342 

"49 

169 

1368 

4828 

9704 

"12150 

9704 

"4828 

1368 

"169 

"7 

42 

"98 

98 

0 

"98 

98 

"42 

7 

39 

"312 

1092 

2184 

2730 

2184 

1092 

"312 

39 

TABLE  1 1. 3  (continued) 
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APPENDIX  III 


VALUES  AND  DESK-COEFFICIENTS 

In  this  appendix  three  tables  are  included.  In  the  first  table 

Table  III.l  the  various  continuous  piecewise  interpolation  formulae  are 

summarized.  In  the  second  table  Table  III. 2  the  values  of  3  =  min  e. . 

1 3 

(3.83)  and  | |ot  F-G |[^  are  given  for  each  formula.  The  value  of  6  is  ob¬ 
tained  using  the  function  VALABETA  while  the  value  of  1 1 a  F-G 1 1 ^  is  calcu¬ 
lated  by  EANORMACPIF.  The  third  table  Table  III. 3  is  the  desk-coefficients 
for  Formula  2.5.2  in  increments  of  .01.  These  coefficients  are  generated 
by  using  the  functions  DESKACOEFF,  POLY,  EFT,  CPI  and  INVERSE.  There  are 
six  coefficients,  Ci,  i  =  1,  2,  ...,6,  for  each  value  of  x  =  ~.5  to  x  =  .5. 
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TABLE  III.l 


SUMMARY  OF  THE  FORMULAE 


FORMULA 

TOTAL  NUMBER 

OF  POINTS  USED 

HIGHEST  ORDER 
OF  CONTINUOUS 
DERIVATIVE 

ORDER  OF  THE 
MULTISTEP 

LOWEST  DEGREE 
OF  POLYNOMIAL 
USED 

2.3.1 

4 

1 

1 

2 

2.5.1 

6 

1 

1 

3 

2.7.1 

8 

1 

1 

3 

2.9.1 

10 

1 

1 

3 

2.11.1 

12 

1 

1 

3 

3.3.1 

5 

1 

2 

2 

3.5.1 

7 

1 

2 

4 

3.7.1 

9 

1 

2 

4 

3.9.1 

11 

1 

2 

4 

3.11.1 

13 

1 

2 

4 

4.3.1 

6 

1 

3 

2 

4.5.1 

8 

1 

3 

4 

5.3.1 

7 

1 

4 

2 

5.5.1 

9 

1 

4 

4 

2.5.2 

6 

2 

1 

4 

2.7.2 

8 

2 

1 

5 

3.5.2 

7 

2 

2 

4 

3.7.2 

9 

2 

2 

6 

3.9.2 

11 

2 

2 

6 

3.11.2 

13 

2 

2 

6 
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TABLE  III  .1  (continued) 


FORMULA 

TOTAL  NUMBER 

OF  POINTS  USED 

HIGHEST  ORDER 
OF  CONTINUOUS 
DERIVATIVE 

ORDER  OF  THE 
MULTISTEP 

LOWEST  DEGREE 
OF  POLYNOMIAL 
USED 

4.5.2 

8 

2 

3 

4 

4.7.2 

10 

2 

3 

6 

5.5.2 

9 

2 

4 

4 

2.7.3 

8 

3 

1 

6 

3.7.3 

9 

3 

2 

6 

3.9.3 

11 

3 

2 

8 
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TABLE  III. 2 

VALUES  OF  3  AND  ||aF  -  G  ||E 


FORMULA 

3  =  min  e .  . 

'  J 

1 1  a F  -  G||e 

2.3.1 

0 

0 

2.5.1 

~  5 . 68434 1 886E- 1 4 

6.200732257E" 

2.7.1 

~2 .8421 70943E~1 3 

4.62490556E‘l 

2.9.1 

1 . 455191 523E~1 1 

2 . 1 0501 9221 E~ 

2.11.1 

"7 . 27595761 4E~1 2 

1.802631396E" 

3.3.1 

0 

0 

3.5.1 

3 . 55271 3679E"l 5 

1 . 26233871 4E” 

3.7.1 

4.618527782E"14 

1 . 1 5774491 4E” 

3.9.1 

~5 .684341 886E- 1 3 

1 . 394752431 E" 

3.11.1 

1 .81 8989404E- 1 2 

5 . 72096741 8E~ 

4.3.1 

"1.534772309E"12 

3.228271377E" 

4.5.1 

~3 . 637978807E- 1 2 

7 . 53061 5877E" 

5.3.1 

"1 . 47792889E” 1 2 

1 . 51 7869634E” 

5.5.1 

~3 .01 2701 2E”l 2 

3.085634286E" 

2.5.2 

”2 . 7284841 05E"1 2 

5.074619662E” 

2.7.2 

"4.365574569E”1 1 

7.939614894E” 

3.5.2 

2 . 273736754E"! 3 

5.0396641 76E" 

3.7.2 

~8. 18545231 6E~1 2 

1 . 759273096E” 

3.9.2 

"2 . 764863893E- 1 0 

5.3941 161 13E" 

3.11.2 

~1 . 1 20497473E- 9 

2 . 531 980834E" 

14 

3 

11 

11 

14 

13 

12 

12 

12 

12 

12 

12 

12 

11 

13 

11 

10 

9 


TABLE  III. 2  (continued) 


FORMULA 

4.5.2 

4.7.2 

5.5.2 

2.7.3 

3.7.3 

3.9.3 


3  =  min  e.. 

LI 

2 .000888344E”1 1 
3.25962901 1E~9 
8 . 622009773E” 1 0 

7.217749953E"9 

8.731149137E'll 

8.003553376E~10 


IhF-  GHe 

5 . 48686651 4E~ 1 1 
6.949643684E"9 
1 .9080448E”9 

1 . 7  89202523E- 8 

1.996080858E"10 
2.201 444096E- 9 


TABLE  III. 3 
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X 

5 . 0E”01 
4  .  9E~01 
4  .  8E~01 
4 . 7E_01 
4  .  6E  01 
4 . 5E  01 
4 . 45’”01 
4 . 3E~01 

4  .  25’_01 
4 . 1E"01 
4 . 0£“01 
3 . 9E~01 
3  .  8E~01 
3  .  72?*" 0 1 
3 . 6Z?”0 1 
3 . SE~ 01 
3  .  42?  01 
3  .  32?”  0  1 
3 . 22?”  0  1 
3 . IE  01 
3 . 0E~01 
2 . 9E~01 
2 . 8£_01 
2 . IE -0 1 
2.6  2?  01 

2  .  52?**0  1 
2.4  E~01 
2 . 3E“01 
2 . 25’”01 
2 . 1E”01 
2  .  02?_ 0 1 
1 . 9E  01 
1 . 82?“01 
1  .  72?_01 
1 . 6£“01 
1 . 5E“01 
1.4E  01 

1  .  32?” 0  1 
1 . 2E“01 
1 . 1E_01 
1.0 E  01 
9  .  02?~0  2 
8 . OE  0  2 
7 . 02?_02 
6 . OE  _0  2 

5  .  02?  02 
4 . 02?”  0  2 

3  .  02?~0  2 

2  .  02?_0  2 
1 . 0E  02 


TABLES  OF  COEFFICIENTS  FOR  FORMULA  2.5.2 


Cl 

1  .  5  612  5  1  12?”l  7 
8  .  2  8  7  9  7  0  62?”  0  4 

1 . 6  4  7  0  8  6  OE ^03 
2 . 4528087E203 
3  .  2440  3  2  02?”03 
4. 018  9  4  5  32?  03 

4 . 7  7  5  8  5  8  0E~ 0 3 

5  .  51319  6  92?~0  3 
6 .  2  2  9  5  0  4  0E~ 0 3 

6  .  9  2  3433  62?~0  3 

7  .  5937500E”03 

8  .  2  3  9  3  2  5  22?”0  3 
8  .  859136  OE  03 
9 . 45226182?“03 
1 .  00178825’”02 
1  .  05552732?~02 
1  .  106  38082?_02 
1  .  154295  02?  ”02 
1  .  19922545’_02 
1  .  2  4113  6  2E~ 0 2 
1  .  2  8  0  0  0  0  02?~0  2 
1  .  31579785’“02 
1  .  3485186E”02 
1  .  3  7  815  9  02?”0  2 
1 . 4047232E_02 
1 . 42822272?_02 
1 . 44867585’_02 

1 . 4  6  61  0  7  82?  ”  0  2 

1 . 4  8  0  5  5  0  42?  _  0  2 
1.4920415E  02 
1  .  5  0  0  6  2  5  02?”  0  2 
1  .  5063506E_02 
1  .  50927365’_02 
1  .  50945432?_02 
1  .  5  0  6  9  5  8  2E  02 
1  .  5  01  8  5  5  52?  ”  0  2 
1 . 494220  82?”02 
1  .  4  8  4 1  3  3  12?”  0  2 
1 . 47167545’_02 
1 . 45693435’_02 
1.440000  OE  02 

1 . 4  2  0  9  6  6  02?~0  2 
1  .  3  9  9  9  2  8  62?”  0  2 
1  .  3  7  6  9  8  7  12?^  0  2 
1  .  3  5  2  2  4  3  22?_  0  2 
1  .  32580082?  02 
1  .  297765  82? ”02 
1 . 2  6  8  2  4  5  9Z?~0  2 
1  .  2  3  7  3  5  0  42?_0  2 
1.2051896E  02 


C  2 

1  .  52655672?”l6 
6  .  5984016£’”03 
1 .  3  0  5  4  0  9  OE” 0  2 

1  .  9  3582602?”02 

2  .  5  5  0  2  7  2  02?~02 

3  .  14  7  9  8  8  3E~  0 2 
3  .  7  2  8  2  7  5  02?~0  2 
4.  29049  01£’”02 
4.83404802?  02 
5  .  3  5  8  418  42?~0  2 

5  .  863125  OS’  02 

6  .  3  4  7  7  4  4  22?~0  2 

6  .  8119040E”02 
7.2552825 E  02 

7  .  6  7  7  6  0  7  0E~ 0  2 
8 . 0786523E”02 
8 . 4582400Z?”02 

8  .  8162367E_02 
9 . 15  2  5  5  3  OE  02 
9 . 46714252?”02 

9  .  7  6  0  0  0  0  0E~ 0 2 
1  .  0031161E”01 
1.0280699E  01 
1  .  0508727£’"01 
1  .  0  715  3  9  2Z?”0  1 
1  .  0900879E”01 
1  .  10  6  5  4  0  5£’”0  1 
1  .  12092212?~01 
1  .  1  3326  082?  01 
1  .  14358792?~01 
1.1519  3  7  52?  01 
1  .  15  8346  5E"01 
1  .  16  2  8  5  4  4E”0 1 
1  .  165  5033E_01 
1  .  16633772?  01 
1  .  1654043Z?~01 
1  .  16  2  7  5  2  0E”0 1 
1  .  15843172?  01 
1  .  152496  32?~01 
1  .  1450003E_01 
1  .  13  6  0  0  0  02?^  0  1 
1  .  12  5  5  5  3  IE _0 1 
1  .  11371892?  01 
1  .  10055772?”01 
1 . 0  861  31  2£’_01 
1  .  0  7  0  5  0  2  02?_  0  1 
1  .  0537  3  352?_01 
1.0  3589012?  01 
1 . 01703682?~01 
9.9723897E  02 


C3 

1  . OOOOOOOEOO 
9.998721 4E  01 
9 . 99477502?”01 
9.988004  52?  01 
9 . 97826  562?”01 
9 . 965425  82?_01 
9  .  9  49  364  22?  01 
9 . 9299713E"01 
9.907148  82?  01 
9  .  8  8080902?”01 
9 . 8508750Z?”01 
9  .  8  17  28022?”01 
9 . 77996802?"01 

9 . 7  3889  18E”01 
9 . 69401462?~01 
9 . 64530862?”01 
9  .  5  9  2  7  5  5  2E”0 1 
9  .  5  36  34472?~01 
9 . 47607582?”01 
9 . 4119558S'”01 
9  .  3  44000  OS’  01 
9  .  272231  5E"01 
9  .  196681  02?“01 
9  .  11738662?"01 
9  .  03439362?”01 
8 . 94775  39E”01 
8  .  8  575  262E”01 
8 . 76  3775  52?“01 
8.66657282?  01 
8  .  56599  512?”01 
8  .  4  6  21  2  5  02?”  0 1 
8  .  35505032?  01 
8 . 24486402?”01 
8  .  1316  6  4  02?_  0  1 
8.0  15552  OE  01 
7  .  8  966  3672?”01 
7 . 7750272E_01 
7  .  65083882?  01 
7  .  5241898Z?”01 
7  .  395  20192?~01 
7.2640000E  01 
7  .  13071162?~01 
6.995467  OE  01 
6. 85839882?”01 
6.719641652  01 
6  .  57933205’”01 

6 . 4  3  7  6  0  8  22?”0  1 
6.29460962?  01 
6  .  1504768E’~01 
6 . 00535135’”01 
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TABLE  III. 3  (continued) 

C  4  c  5  c  6 


5  . 

,  05 

01 

1  . 

4  . 

,  95" 

"01 

6  . 

4. 

,  85" 

"01 

1  . 

4  . 

,  75" 

"01 

2  . 

4  . 

,  65" 

"01 

2  . 

4. 

,  55" 

"01 

3  . 

4. 

,  45" 

"01 

4  . 

4. 

,  35" 

"01 

5  . 

4. 

,  25" 

"01 

5  . 

4  . 

,  15" 

"01 

6  . 

4. 

,  05" 

"01 

7  . 

3  . 

,  95" 

"01 

8  . 

3  . 

,  85" 

"01 

9  . 

3  . 

,  75" 

"01 

1  . 

3  . 

,  65 

"01 

1  . 

3  . 

,  5 5" 

"01 

1  . 

3. 

,  45" 

"01 

1  . 

3. 

,  3 5" 

"01 

1  . 

3. 

,  25" 

"01 

1  . 

3  . 

,  15" 

"01 

1  . 

3  . 

.  05" 

"01 

1  . 

2  . 

,  95" 

"01 

1  . 

2  . 

,  85" 

"01 

1  . 

2  . 

,  75" 

"01 

2  . 

2  . 

,  65" 

"01 

2  . 

2  . 

,  55" 

"01 

2  . 

2  . 

,  45" 

"01 

2  . 

2. 

,  35" 

"01 

O 

^  . 

2  . 

,  25" 

"01 

2  . 

2  . 

,  15" 

"01 

2  . 

2  . 

,  05" 

"01 

2  . 

1  . 

,  95" 

"01 

3  . 

1  . 

.  85" 

"01 

3  . 

1  . 

.  75" 

"01 

3  . 

1  . 

,  65" 

*01 

3  . 

1  . 

.  55" 

"01 

3  . 

1  . 

,  45" 

"01 

3  . 

1  . 

,  35" 

"01 

3  . 

1  . 

,  25" 

"01 

4  . 

1  . 

,  15" 

"01 

4  . 

1  . 

.  05" 

"01 

4. 

9  . 

.  05" 

"02 

4. 

8  . 

.  05" 

"02 

4. 

7  . 

.  05" 

"02 

4  . 

6. 

.  05" 

"02 

4. 

5  . 

,  05" 

‘02 

5  . 

4  . 

.  05" 

"02 

5  . 

3  . 

.  05" 

"02 

5  . 

2  . 

,  05" 

"02 

5  . 

1  . 

.  05" 

‘02 

5  . 

0  9  6  3  4  5  25^1 5 

7  3  6  0  31  95*0  3 
3621  18  05*0  2 
06701165*02 
78963205*02 
5  31210  95*0 2 

2  9  2  8  6  6  05*0 2 
0  7  5  6  0  3  35*0  2 

8  8  0  3  2  0  05*0 2 
70780675 ”02 
55875005_02 

4  3  3  7  3  5  15*0  2 

3  3  3  2  4  8  05*0  2 
02576785'  01 
12073225*01 
21823835_01 
31829765^01 
42091305  01 
52607905*01 
63378195_01 

7  4  4  0  0  0  05  ^  0 1 

8  5  6  7  0  4  25  ” 0 1 
9718578 5  01 
0  8  9  417  05  *  0 1 
20933125  01 
3  315  4  3  05*0 1 
45598865*01 

5  8  2  5  9  8  25*01 
7112  9  6  05^0 1 
84200055*01 

9  7  4  6  2  5  05*0 1 
10907745  01 

2  4  5  2  6  0  85*0 1 
38307375*01 

5  2  2  410  25*0 1 
66316025  01 
80520965*01 
94844095*01 
0  9  2  7  3  3  05  *01 
23796175*01 

3  8  4  0  0  0  05  *  0 1 
530718 15*0 1 

6  7  7  9  8  3  85  _  0 1 
8256629 5  01 
9  7  3  619  25*01 
12171485*01 
26981065_01 
41776615^01 
5  6  5  4  4  0  05” 0 1 
71269045  01 


5 . 551115 15^17 
8 . 38849695*04 
1 .69393005*03 
2.  5  7  2  5  91 65*0  3 
3. 48160005^03 

4 . 4  2  714  8  45_  0  3 
5.41487005  03 

6 . 4  49850  35*0  3 
7 . 53664005*03 
8 . 679  26  7  25*03 

9 . 8  812  5  0  05  ^  0  3 
1 . 11456095*02 
1 . 24748805*02 
1 . 38711 265*02 
1 . 53359505*02 
1  .  6  8  7  0  5  0  85**0  2 
1  .  8  4  7  5  5  2  05*0  2 
2.0151285 5  02 
2 . 1897  69  05*02 
2 . 37142275^02 
2.56000005  02 
2 . 75  537435*02 
2 . 95738305*02 
3  .  16  5  8  2  8  55*0  2 

3 . 3  8  0  4  8  0  05*0  2 

3 . 6  010  74  25*0  2 

3. 8  2  7  317  05*0  2 
4.05888445  02 

4. 2  9  5  4  2  4  05” 0  2 
4.5365561 5  02 

4. 7  818  7  5  05*0  2 
5 . 0  3  0  9  5  0  35*0  2 
5  .  2  8  3  3  2  8  05*0  2 
5 .  5  3  8  5  3  2  05^0  2 
5.79606505  02 
6 . 0  5  5  410  25*0  2 
6.31603205  02 
6 . 57737785*02 

6 . 8  3  8  8  7  9  05*0  2 
7 . 0  9  9  9  5  2  05*0  2 

7 . 3  6  0  0  0  0  05  ^  0  2 
7 . 61841375*02 
7.8745730 5  02 
8 . 12784795*02 
8.37760005  02 
8 . 6  2  318  3  65~0  2 
8 . 86394705^02 
9 . 0  9  9  2  3  3  85~  0  2 

9 . 3  2  8  3  8  4  05~0  2 
9 . 55073555“02 


6.0715322 5  18 
2 . 86687505"07 
2. 254000 05 "06 
7.47506255*06 
1 . 74080005"05 
3 . 33984385"05 
5 . 66820005"05 
8 . 83868135"05 
1 .29536005^04 
1 . 81050195^04 
2 . 43750005"04 
3 . 18358565^04 
4 . 05504005_04 
5.0572194 5  04 
6 . 19  4  5  8  0  05~04 
7.47070315  04 
8 . 88832005"04 
1 . 04493375"03 
1 . 2 1 548605~03 
1.40052215  03 
1 . 60000005"03 
1 . 8  13  8  0  5  45~0  3 
2 . 04175405"03 
2.28359385  03 
2 . 5  3  9  0  0  8  05~  0  3 
2.80761725  03 
3 . 08898205"03 
3.3826056 5  03 
3 . 68793605"03 
4 . 00436895*03 
4. 33125005*03 
4 . 66787735*03 
5 .01350405*03 
5.36734075  03 
5 .72855805*03 
6 .09628915*03 
6 . 46963205*03 
6 . 84765245*03 
7 .22938605*03 
7 . 61384085*03 
8 . 00000005*03 
8 . 38682425*03 
8 .77325405*03 
9.1582126 5  03 
9 . 54060805*03 
9.9193359 5  03 
1 .029328  25^02 
1 . 066 13245_02 
1 . 10  2  2  3  3  65*0  2 
1.1375188 5  02 
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TABLE  III. 3  (continued) 


* 

Cl 

C  2 

C  3 

5  .  2#” 

18 

1  .  1718750#” 

02 

”9.7656250 #” 

02 

T 

5,8593750#" 

1  .  0 #” 

02 

1  .  1375188 E~ 

02 

“9  .  55073552? 

02 

5.7126904# 

2  .  0 E~ 

02 

1 . 1022336 E~ 

02 

"9  .  328  3840 2?” 

02 

5  .  5  6  5  4  4  0  0Z?” 

3. 0 E~ 

02 

1  .  0  6  613  2  4#” 

02 

”9  .  09923382? 

02 

5 .4177661#" 

4 . 0#” 

02 

1  .  0293282 E~ 

02 

8  .  86394702?” 

02 

5  .  2698106#” 

5  .  0  E 

02 

9 . 9193359#“ 

03 

"8 .62318  36 2?” 

02 

5 . 1217148#” 

6 . 0  E 

02 

9  .  5406080#" 

03 

"8  .  37760002?" 

02 

4 . 9736192#” 

7  .  0#” 

02 

9  .  1582126 E~ 

03 

”8  .  1278479 2?" 

02 

4.8256629# 

8  .  0 E~ 

02 

8  .  7732540#” 

03 

”7  .  87457302? 

02 

4.6779838# 

9  .  0 E~ 

02 

8  .  3868242#” 

03 

”7 .61841372?" 

02 

4 . 5307181#” 

1  .  0 #" 

01 

8  .  0000000#" 

03 

”7  .  36000002?” 

02 

4 . 3840000#” 

1 .  1#” 

01 

7 .6138408#" 

03 

”7  .  0  9  9  9  5  2  02?” 

02 

4 . 2379617#” 

1 .  2E~ 

01 

7  .  2  2  9  3  8  6  02?” 

03 

"6  .  8  3  8  8  7  9  02?” 

02 

4 . 0927330#“ 

1  .  3 #” 

01 

6  .  8476524# 

03 

”6  .  5773778#” 

02 

3 . 9484409#” 

1 . 4  E 

01 

6 . 46963202? 

03 

“6  .  316  0  3  2  02?” 

02 

3 .8052096#" 

1  .  5 E~ 

01 

6 . 09628912?" 

03 

"6  .  05541022? 

02 

3 .6631602#” 

1  .  6 E~ 

01 

5  .  7285580 2?" 

03 

"5  .  7  9  6  0  6  5  02?" 

02 

3 . 5224102# 

1  .  IE 

01 

5  .  36734072? 

03 

"5  .  5  3  8  5  3  2  02?” 

02 

3 . 3830737#" 

i .  8#” 

01 

5 . 01  350402?” 

03 

”5  .  283  328027" 

02 

3 .2452608#" 

1  .  9 E~ 

01 

4 . 6678773#” 

03 

”5  .  0  3  0  9  5  0  32?" 

02 

3  .  1090774# 

ro 

• 

o 

tg 

i 

01 

4 . 3  312  5  0  02? 

03 

"4 . 7  818  7  5  OE 

02 

2 . 9746250#" 

2  .  IE 

01 

4  .  0  0  4  3  6  8  92?” 

03 

”4 . 5365561#” 

02 

2  .  8420005# 

2 . 2E~ 

01 

3 . 68793602?" 

03 

“4 . 2  9  5  4  2  4  0E 

02 

2 .7112960#" 

2  .  3  E 

01 

3  .  3826056 E 

03 

”4. 05  8  8  8442?" 

02 

2 .5825982# 

2 . 4# 

01 

3  .  0  8  8  9  8  2  02?” 

03 

”3.  8  2  7  317  02?” 

02 

2 .4559886#” 

2  .  5#” 

01 

2 .8076172 2?" 

03 

”3. 60107422?" 

02 

2 . 3315430#" 

2 . 6 #” 

01 

2  .  5  3  9  0  0  8  02?" 

03 

”3  .  38  0480 02?” 

02 

2  .  2093312# 

2 . 1E~ 

01 

2 .2835938 2?” 

03 

"3  .  16582852?" 

02 

2 . 0894170#” 

2  .  8 E~ 

01 

2  .  041754 02?” 

03 

”2  .  9  5  7  3  8  3  02?” 

02 

1 . 9718578#" 

2  .  9 E~ 

01 

1 .8138  05  42?” 

03 

”2  .  7553743#” 

02 

1 . 8567042#” 

3  .  0 #" 

01 

1 .60000002?” 

03 

”2  .  56000002? 

02 

1 . 7440000#” 

3  .  IE" 

01 

1 .4005221#” 

03 

”2  .  3714  2  272?” 

02 

1 .6337819#" 

3  .  2 E~ 

01 

1  .  21  548602?” 

03 

"2  .18976902? 

02 

1.5260790# 

3  .  3 E~ 

01 

1  .  04493372?" 

03 

"2  .  0151285 2?” 

02 

1 .4209130#” 

3 . 4# 

01 

8  .  8  8  8  3  2  0  02?” 

04 

”1. 8475520 E~ 

02 

1 . 3182976#” 

3  .  5  E 

01 

7 . 4707031#” 

04 

"1.6870508 #” 

02 

1 .2182383#" 

3 . 6#” 

01 

6  .  1  9458  0  02? 

0  4 

1  .  533595  02? 

02 

1 . 1207322#” 

3  .  IE 

01 

5 . 0572194#” 

04 

"1  .  3871  1262?” 

02 

1 .0257678 #" 

3  .  8  E 

01 

4 . 0  5  5  040  02? 

04 

1 .  2474880#” 

02 

9  .  3332480# 

3  .  9 E~ 

01 

3  .  183  58562? 

04 

”1 .  11456092?” 

02 

8.4337351# 

4 . 0 E~ 

01 

2 .437  5  0002?” 

04 

"9  .  8812500 2?” 

03 

7 . 5587500#” 

4 . 1  E~ 

01 

1 . 8105019#” 

04 

8 . 6792672# 

03 

6 . 7078067#" 

4 . 2E~ 

01 

1  .  2  9  5  3  6  0  02?” 

04 

”7  .  5  3  6  6  4  0  02?” 

03 

5 . 8803200#" 

4 . 3E~ 

01 

8  .  8386813#” 

05 

“6 . 4  4  9  8  5  0  3E~ 

03 

5 . 0756033#” 

4 . 4 E~ 

01 

5 .6682000#” 

05 

"5 .41487002?" 

03 

4 .2928660#” 

4 . 5  E 

01 

3  .  3398438#” 

05 

”4 . 42714842? 

03 

3  .  5312109# 

4 . 6  E 

01 

1  .  74080002? 

05 

"3. 48160002? 

03 

2 .7896320#” 

4 . 7#" 

01 

7 .4750625#” 

06 

”2  .  57259162?” 

03 

2 . 0670116#” 

4 . 8  E 

01 

2  .  2  5  4  0  0  0  0E~ 

06 

"l  .  6939300# 

03 

1 . 3621180#" 

4 . 9  E 

01 

2 .8668750 E~ 

07 

8  .  38849692?” 

04 

6 .7360319#” 

5  .  0 E~ 

01 

"3.4694470 E 

18 

2  .  7755  5762? 

17 

1 .4571677#” 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

01 

02 

02 

02 

02 

02 

02 

02 

02 

02 

02 

02 

03 

15 
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X 


C  4 


C  5 


C6 


5  .  2  E 

18 

5  . 

1  .  0E~ 

02 

6  . 

2  .  OE 

02 

6  . 

3  .  0 E~ 

02 

6  . 

4  .  OE 

02 

6  . 

5  .  OE 

02 

6  . 

6  .  OZT 

02 

6  , 

7  .  OE 

02 

6  . 

fci 

o 

• 

CO 

02 

6  . 

9  .  OE 

02 

7  . 

1  .  OE" 

01 

7  . 

1  .  IE" 

01 

7  . 

1  .  2 E~ 

01 

7  . 

1  .  3 E~ 

01 

7  . 

1 . 4Z?" 

01 

7  . 

1  .  5 E~ 

01 

7  . 

1  .  6  E 

01 

8  . 

1  .  IE" 

01 

8  . 

1  .  8 E~ 

01 

8  . 

1  .  9 E~ 

01 

8  . 

2  .  OE" 

01 

8  . 

2 . 1ZT 

01 

8  . 

2  .  2 E~ 

01 

8  . 

2  .  3 E~ 

01 

8  . 

2 . 45’” 

01 

8  . 

2  .  5 E~ 

01 

8  . 

2 . 6  E 

01 

9  . 

2  .  7  E 

01 

9  . 

2  .  8  E 

01 

9. 

2  .  9 E~ 

01 

9. 

OJ 

• 

o 

l 

01 

9  . 

3  .  IE 

01 

9  . 

3 . 2 E~ 

01 

9  . 

3  .  3 E~ 

01 

9. 

3 . 4Z? 

01 

9  . 

3  .  5  E 

01 

9  . 

3 . 6 E~ 

01 

9  . 

3  .  IE 

01 

9  . 

3  .  BE" 

01 

9  . 

3  .  9 E~ 

01 

9  . 

4  .  OE 

01 

9  . 

4  .  IE" 

01 

9  . 

4 . 2 E~ 

01 

9. 

4  .  BE 

01 

9  . 

4 . 4£Z 

01 

9. 

4 . 5  E 

01 

9. 

4 . 6  E 

01 

9  . 

4  .  IE" 

01 

9  . 

4  .  BE" 

01 

9  . 

4 . 9 E~ 

01 

9  . 

5  .  0E~ 

01 

1  . 

8  5  9  3  7  5  OS’  01 

00535135-01 

15  047  6  8Z?"~  0 1 
2  94609  SE~ 0 1 
4376082E"01 
57933205"01 
719641 6Z?~  0 1 
8  5  8  3  9  8  8Z7”  0 1 
9954670Z?  01 
13  0  711 6E~  0 1 
2  6  4  0  0  0  OE -0 1 
3952019Z?  01 
5241898£’~01 
6508388£~01 

7  7  5  0  2  7  2Z?”  0 1 

8  9  6  6  3  6  7Z?  01 
015  5  5  2  6E~ 0 1 

1  3  1 6  6  4  OZ?”  0 1 

2  4  4  8  6  4  0E~  0 1 

3  5  5  0  5  0  3E  01 
46  2  125  OZ?" 01 
5659951£~01 
6665728£'"01 

7  6  3  7  7  5  5Z7~  0 1 

8  5  7  5  2  6  2Z7~  0 1 
94775395,_01 
03439365~01 
1 1  7  3  8  6  6Z7~  0 1 

1  96681 OE  _ 0 1 

2  7  2  2  31  5  £7  ~  0  1 

3  44000  OE  01 
411  9558£’~01 
4760758£’~01 
5363447Z?"01 
5  9  2  7  5  5  2Z7~  0  1 
64530865’_01 
6940146£’_01 
7388918S’_01 
77996805_01 
81728025"01 
8  5  0  8  7  5  0E~  0 1 

8  8  0  8  0  9  0E~ 0 1 
90714885’~01 
9299713£’_01 

9  4  9  3  6  4  2Z7  01 
9  6  5  4  2  5  SE~ 0 1 
9  7  8  2  6  5  BE^O 1 
98800455’_01 
9  9  4  7  7  5  OE  01 
9  9  8  7  21 4Z?*"  0 1 
0000000500 


9  .  7  6  5  6  2  5 OE  02 
9.  97  238  97Z7~02 
1  .  0170368Z?"01 
1 . 0  3  5  8  9  01E~  0 1 
1  .  0  5  3  7  3  3  SE~  0 1 
1 . 070502  0E~  0 1 
1  .  08613125~01 
1 . 1005577Z?“01 
1  .  1137189£’"01 
1  .  1 2  5  5  5  3 1Z?”  0  1 
1  .  13600005~01 
1  .  14  5  0  0  0  3E~ 0 1 
1 . 15249635~01 
1  .  15843175~01 
1  .  16  2  7  5  2  OS’”  0 1 
1.16540435  01 
1  .  16633775_01 
1  .  16  5  5  0  3  35”  0  1 
1 .  16  2  8  5  4  4#  01 
1  .  1583465£’201 
1.  15193755  01 
1 . 14  3  5  8  7  9  5~ 0 1 
1  .  13326085~01 
1.12092215  01 
1  .  10654055~01 
1  .  0  9  0  0  8  7  95” 0 1 
1. 071  5  39  2E~ 0 1 
1 . 05087275“01 

1  .  0  2  8  0  6  9  95”  0  1 
1.00311615  01 
9  .  7  6  0  0  0  0  05”  0  2 
9. 4671425Z?_02 
9  .  15  2  5  5  3  OE  02 
8  .  8162367£’~02 
8 . 4582400£^02 
8.078652  3E  02 
7 . 6776070£"02 
7.2552825E1  02 
6  .  8119  0  4  027~  0  2 
6  .  3477442£’“02 
5  .  8631250£’~02 
5  .  3  5  8  418  4Z?”  0  2 
4 . 8  3  4  0  4  8  0E~  0 2 
4 . 2  9  049  0  IE1”  0  2 
3  .  7282750Z?"02 
3  .  14  7  9  8  8  3Z?~  0  2 

2  .  5  5  0  2  7  2  0E~  0 2 
1  .  9  3  5  8  2  6  0E~  0 2 
1  .  3  0  5  4  0  9  OE1"  0  2 
6 . 598401 6E~ 0 3 
4 . 4408921£"16 


1  .  17187505’  02 
1  .  2051896E’~02 
1 . 2373504£’"02 
1 . 26824595’_02 
1 . 29776585’_02 
1.  325800851  02 
1  .  3  5  2  2  4  3  2E~  0 2 
1  .  376987 1E~  0 2 
1  .  3  9  9  9  2  8  6E~ 0 2 
1 . 4209660S’“02 

1 . 4  4  0  0  0  0  0E~  0 2 
1 . 4569343£’”02 
1 . 471675 45’"' 0 2 
1 . 48413315’"02 
1 . 494220827  02 
1  .  5018  55  5£’“02 
1  .  50695825’~02 
1  .  5  0  9  4  5  4  35’’*  02 
1  .  50927  36£’_02 
1  .  5  0  6  3  5  0  627** 0 2 
1  .  5  0062  505’_02 
1 . 492  0415£’”02 
1 . 4805504£’"02 
1 . 466  10785’“02 
1 . 44867585’"02 
1 . 42822275"02 
1 . 4047232E’”02 
1 . 378159  0E~ 02 
1  .  3485  186Z?"02 
1  .  31579785"02 
1  .  2  8  0  0  0  0  0E~ 0 2 
1 . 24113625-02 
1  .  19  9  2  2  5  4Z?~  0  2 
1  .  15  4  2  9  5  05” 0  2 
1  .  10  6  3  8  0  QE~  0 2 
1  .  0555273E’~02 
1  .  0017882£’“02 
9 . 45226185’"03 
8 .  ssgiseofi^os 
8  .  23932525’~03 

7 . 5  9  3  7  5  0  0E~  0 3 
6  .  9  2  3  4  3  3  0E~ 0 3 
6  .  2  2  9  5  04  0E~  0 3 
5  .  51319697T03 
4.77585  8  OE  03 
4 . 01894535’"03 
3  .  2  4  4  0  3  2  05” 0  3 
2 . 45280875'"03 
1 . 64708605’"03 
8 . 287  97  065~04 
1  .  9  9  4  9  3  2  OZ?”  1 7 


APPENDIX  IV 


NUMERICAL  RESULTS  OF  THE  TEST  FUNCTIONS  AND  PROCEDURE 

In  this  appendix  three  tables  are  given.  Table  IV. 1  gives  the 
number  of  multiplications  required  in  the  application  of  each  coefficient 
matrix  F  in  (3.53)  by  the  procedure  given  in  section  5.2  of  Chapter  V. 
The  number  of  multiplications  required  for  the  set  of  coefficients  of  H 
corresponds  to  that  required  in  calculating  all  the  coefficients  e^  of 
(5.4)  for  each  coefficient  matrix  F  of  (3.53).  The  number  of  multi¬ 
plications  in  one  evaluation  of  H  is  that  required  in  evaluating  a 
polynomial  of  the  form  (5.6).  In  Tables  IV. 2  and  IV. 3  are  given  the 
maximum  relative  errors  for  (5.1)  and  (5.2)  respectively .  The  two  last 
tables  were  obtained  using  the  functions  APPLICACPIF,  POLY,  CPI,  and 
INVERSE.  The  results  in  Table  IV. 2  are  for  a  constant  interval  i.e. 
[1.025,2],  while  those  in  Table  IV. 3  are  for  the  intervals  [10.025,11] 
and  [25.025,26]. 
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TABLE  IV. 1 


NUMBER  OF  MULTIPLICATIONS  REQUIRED 


FORMULA 

THE  SET  OF  COEFFICIENTS  OF  H 

ONE  EVALUATION 

2.3.1 

8 

3 

2.5.1 

12 

3 

2.7.1 

16 

3 

2.9.1 

20 

3 

2.11.1 

24 

3 

3.3.1 

11 

4 

3.5.1 

15 

4 

3.7.1 

19 

4 

3.9.1 

23 

4 

3.11.1 

27 

4 

4.3.1 

18 

5 

4.5.1 

24 

5 

5.3.1 

22 

6 

5.5.1 

28 

6 

2.5.2 

18 

5 

2.7.2 

24 

5 

3.5.2 

22 

6 

3.7.2 

28 

6 

3.9.2 

34 

6 

3.11.2 

40 

6 

4.5.2 

32 

7 

4.7.2 

40 

7 
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TABLE  IV. 1  (continued) 

FORMULA 

THE  SET  OF  COEFFICIENTS  OF  H  ONE  EVALUATION  OF  H 

5.5.2 

37  8 

2.7.3 

32  7 

3.7.3 

37  8 

3.9.3 

45  8 

TABLE  IV. 2 


THE  LARGEST  RELATIVE  ERROR  FOR  x\  i  =  6,7,10 


FORMULA 

2.3.1 

2.5.1 

2.7.1 

2.9.1 

2.11.1 

3.3.1 

3.5.1 

3.7.1 

3.9.1 

3.11.1 

4.3.1 

4.5.1 

5.3.1 

5.5.1 

2.5.2 

2.7.2 

3.5.2 

3.7.2 

3.9.2 

3.11.2 

4.5.2 

4.7.2 


6 

7 

10 

X 

X 

X 

4 . 96431 E~3 

1 . 03283E” 2 

4. 56391 E"2 

3.29647E"4 

7.59429E”4 

-3 . 36421 E” 3 

2.96472E"4 

6.92587E”4 

4. 34825E”3 

2.96472E"4 

6 . 92  587E- 4 

4. 1 83  79E” 3 

2.96472E”4 

6.92587E"4 

4. 1 8507E” 3 

9 . 1 2287E~3 

1 .65640E~2 

6 .601 24E”2 

"1 . 20968E”4 

~  4 . 4001 7E~  4 

”6 . 23479E" 3 

7 . 28776E”5 

2.27503E"4 

2.’93142E"3 

7. 28776E”5 

2 . 22446E”4 

2 . 27444E”3 

7.28776E'5 

2.22446E”4 

2 . 28882E"3 

1 .080&9E”2 

2.00979E”2 

9 . 59634E” 2 

~2.42844E”4 

"8.99391E”4 

"1 . 30269E"2 

6.85474E'3 

1 . 2461 1 E”2 

4. 89673E"2 

*1.31 759E”4 

”4.691 75E”4 

”6.28779E”3 

'1 .2301 7E”4 

"5.00746E”4 

” 7 . 54654E"3 

”8.7791 5E”8 

"4.90022E”6 

4 . 65282E” 4 

”3.00247E”4 

"1.07975E"3 

"1 . 33627E” 2 

2 . 5847 1 E~ 1 5 

”2.10842E”5 

9.011 92E”4 

2 . 71 395E”1 5 

”1.95108E”7 

"4. 4351 4E"5 

~  2 . 94344E” 1 5 

”1 . 951 08E” 7 

” 2 . 03496E " 5 

"4.79520E”4 

”1 . 82021 E"3 

” 2 . 74970E ” 2 

"2 . 721 1 6E~1 5 

1 . 50320E”5 

2.33937E"3 

TABLE  IV. 2  (continued) 


FORMULA  £ 

5.5.2  ”2 .52548E”4 

2.7.3  "l . 44329E- 1 5 

3.7.3  4 . 45544E- 1 5 

3.9.3  3 . 70028E”! 5 


7 

10 

X 

X 

”9.161 53E"4 

”1 . 291 37E~2 

"5.50463E”6 

5 . 56570E"4 

"2.37197E"5 

1.17684E”3 

3.58047E  15 

2.53198E  5 

Note:  Interval  used  was  [1.025,2]  in  increments  of  .025 


TABLE  IV. 3 


THE  LARGEST  RELATIVE  ERROR  FOR  ex 


FORMULA 

INTERVAL  [10.025,11] 

INTERVAL  [25.025,26] 

2.  3  .1 

"2 . 22957E” 7 

”1 . 00527E”8 

2.5.1 

"2.46006E'9 

”4. 73603E”! 1 

2.  7  .1 

"2.50947E”9 

”4 .75101 E'l 1 

2.  9  .1 

”2 . 50920E- 9 

"4.75101 E” 1 1 

2.11  .1 

"2.50920E"9 

”4.75101 E” 1 1 

3.  3  .1 

"4.20991 E”7 

"1 .9631 6E"8 

3.  5  . 1 

4 . 47  228E- 1 0 

3 . 28966E” 1 2 

3.  7  .1 

"2 . 1 4235E"1 0 

"1 . 67370E” 1 2 

3.  9  .  1 

”2 . 1 21 27E"1 0 

”1 .671 09E” 1 2 

3.  1 1  .1 

”2.1 21 43E”1 0 

"1 . 671 02E" 1 2 

4.  3  .  1 

"2.96954E"7 

”1 . 34385E"8 

4.  5  .  1 

4.40824E"l 0 

3 . 2 1 082E” 1 2 

5.  3  .  1 

"3 . 29027E”7 

"1 . 531 83E” 8 

5.  5  .  1 

4.88252E"10 

3.65441E"12 

2.  5  .  2 

4.88991 E”10 

3 . 551 96E” 1 2 

2.  7  .  2 

"1 .91871 E'l 2 

2 . 81 353E" 1 5 

3.  5.2 

9 . 47390E- 1 0 

7 . 1 3548E” 1 2 

3.  7  .  2 

”3 . 69734E"1 2 

5 . 35420E” 1 5 

3.  9  .  2 

1 .04149E"13 

1 . 64445E" 1 5 

3.  11  .2 

7.71 296E- 1 4 

1 . 58609E”! 5 

4.  5  .  2 

7 . 901 88E”1 0 

5.7211 OE” 1 2 

C\J 

r-^ 

”3.211 32E"1 2 

4.1 061 8E" 1 5 

5.  5  .  2 

8 . 64365E” 1 0 

6 . 4391 4E” 1 2 
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TABLE  IV. 3  (continued) 

FORMULA 

INTERVAL  [10.025,11]  INTERVAL  [25.025,26] 

2,  7.3 

"2 . 391 1 0E"1 2  ”3 . 1 01 31 E~1 5 

CO 

e'¬ 

en 

~4.21270E~1 2  5.97381 E~1 5 

3.  9  o  3 

3 . 34744E”! 4  3.47573E'15 

Note:  Increments  used  were  .025 


